Transport_Interface.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - Transport_Interface
%    Copyright (C) 2009-2015 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 Transport_Interface.m
%> @brief A class to evaluate the Dyson equation and to calculate the scattering matrix for equilibrium processes.
%> @} 
%> @brief A class to evaluate the Dyson equation and to calculate the scattering matrix for equilibrium processes.
%%
classdef Transport_Interface < handle & Messages


    properties ( Access = protected )
        %> The energy to be used in the calculations.
        E
        %> An instance of the structure param
        param
        %> The calculated Greens function.
        G
        %> an instance of #junction_sites describing the sites in the calculated Green operator
        junction_sites
        %> The calculated scattering matrix.
        Smatrix
        %> Numbers of the open channels in the leads.
        nyitott_csatornak
        %> Array ofstructures open_channels.
        open_channels        
        %> The matrix of the calculated conductance.
        Conductance_Matrix
        %> A function handle of the Dyson equation.
        Dysonfunc
        %> An instance of class #CreateHamiltonians or its subclass
        CreateH
        %> An instance of class #Peierls.
        PeierlsTransform
        %> An array of classes #Lead (or its subclasses)
        Leads
        %> list of optional parameters (see http://www.mathworks.com/help/matlab/ref/varargin.html for details)
        varargin
    end




methods

%% Contructor of the class
%> @brief Constructor of the class.
%> @param Opt An instance of the structure #Opt.
%> @param param An instance of the structure #param.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
function obj = Transport_Interface(E, Opt, param, varargin)
    obj = obj@Messages( Opt );
    
    obj.E = E;
    obj.param = param;    
    obj.varargin = varargin;
    
    obj.Initialize();

    
end

%% ScatterCalc
%> @brief Calculates the effective (decimated) Hamiltonian of the scattering region. (Obsolete)
function ScatterCalc( obj )
        
    obj.display('Creating Hamiltonian for Scattering region')
    if isempty(obj.CreateH)
        obj.CreateH = CreateHamiltonians(obj.Opt, obj.param);
    end
    obj.CreateH.CreateScatterH();

    if obj.Opt.magnetic_field == 1
        obj.display('Applying magnetic field in scattering region')
        if isempty( obj.PeierlsTransform )
            obj.PeierlsTransform = Peierls(obj.Opt, 'param', obj.param);
        end
        obj.PeierlsTransform.PeierlsTransform(obj.CreateH, 'Hscatter');
    else
        obj.PeierlsTransform = [];
    end

    if ~isempty(obj.param.scatter.Overlap_in_Scatter) && obj.param.scatter.Overlap_in_Scatter && ~(isempty(obj.Opt.type_of_scanning) ) && ~obj.Opt.Just_Create_Hamiltonians
        obj.display('Applying overlap integrals in scattering region')
        obj.ApplyOverlapInScatter( obj.CreateH, obj.E )
    end

    % Deciamtion
    if obj.Opt.Decimation >= 1  && ~obj.Opt.Decimate_only_Dyson && ~obj.Opt.Just_Create_Hamiltonians
        ido = cputime();
        obj.display('Applying decimation procedure in scattering region')
        Decimation_Procedure = Decimation(obj.Opt); 
        Decimation_Procedure.DecimationFunc(obj.E, obj.CreateH, 'Hscatter', 'kulso_szabfokok');
        obj.display(['A decimalasnal eltelt ido:', num2str(cputime()-ido)]);
    else
        obj.display('No decimation calculation is performed in the scattering region')
    end

end

%% LeadCalc
%> @brief Invokes the function #SurfaceGreenFunctionCalculator for each lead.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'shiftLeads' A real number. If given, the on-site potentials of the sites in the lead are shifted by this value.
%> @param 'transversepotential' A function handle pot = f( #coordinates ) or pot=f( #CreateLeadHamiltonians, Energy) of the transverse potential applied in the lead. (Instead of #CreateLeadHamiltonians can be used its derived class)
%> @param 'coordinates_shift' An integer. If given, the coordinates of the sites in the lead are shifted by coordinates_shift*lattice vector.
%> @param 'gauge_field' Function handle S = f( x,y) of the gauge transformation on the Hamiltonians of the leads. (S is a N x 1 vector, where N is the number of the points given by the x and y coordinates.)
%> @param 'createCore' Set 1 for creating the class instances #Lead but omitting any further calculations, or false (default) to continue with the calculations.
%> @param 'SelfEnergy' Logical value. Set true to calculate the self energies of the leads. (default is false)
%> @param 'SurfaceGreensFunction' Logical value. Set true to calculate the surface Green functions of the leads. (default is true)
%> @param 'leads' Array of lead identification numbers for which the calculations should be performed. (Default is a list of all leads)
%> @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 'CustomHamiltonian' An instance of class #Custom_Hamiltonians or its derived class to create custom Hamiltonians in the leads.
%> @param 'Just_Create_Hamiltonians' Logical value. Set true to create the Hamiltonians of the leads, but stop further calculations. Default value is false.
%> @param 'q' The tranverse momentum for transverse computations.
%> @return Returns with an array of the created #Lead instances.
function lead_ret = LeadCalc( obj, varargin )
        
    p = inputParser;
    p.addParameter('shiftLeads', zeros(length(obj.param.Leads),1));
    p.addParameter('transversepotential', []); 
    p.addParameter('coordinates_shift', zeros(length(obj.param.Leads),1) );
    p.addParameter('gauge_field', [] );
    p.addParameter('createCore', 0);
    p.addParameter('SelfEnergy',false); % set true to calculate the self energy of the semi-infinite lead
    p.addParameter('SurfaceGreensFunction', true );% set true to calculate the surface Greens function of the semi-infinite lead
     p.addParameter('leads', 1:length(obj.param.Leads) );% list of leads to be calculated
    p.addParameter('leadmodel', []); %function handle for an individual physical model for the contacts
    p.addParameter('CustomHamiltonian', []);  
    p.addParameter('Just_Create_Hamiltonians', 0);
    p.addParameter('q', [] ) %transverse momentum
    p.parse(varargin{:});
    
    shiftLeads             = p.Results.shiftLeads;
    transversepotential              = p.Results.transversepotential;
    coordinates_shift      = p.Results.coordinates_shift; % shifts coordinates by "coordinates_shift" times of a lattice vector.
    gauge_field            = p.Results.gauge_field; % if a valid function handle is given, it is used to perform gauge transformation on Lead Hamiltonians
    createCore             = p.Results.createCore;
    SelfEnergy             = p.Results.SelfEnergy;
    SurfaceGreensFunction  = p.Results.SurfaceGreensFunction;
     leads                  = p.Results.leads;
    Just_Create_Hamiltonians = p.Results.Just_Create_Hamiltonians;
    q                      = p.Results.q;
    CustomHamiltonian      = p.Results.CustomHamiltonian;
    leadmodel              = p.Results.leadmodel;
        
     if isempty(obj.Leads)
         obj.Leads = cell(length(obj.param.Leads),1);
     end
    obj.display(' ')
    obj.display('Creating surface_Greens_function interfaces:')
    obj.display(' ')
    idx = 1;
    while idx <= length(leads)
        obj.display('------------------------------')
        obj.display(['Lead: ',num2str(leads(idx))]) 
        obj.Leads{leads(idx)} = obj.SurfaceGreenFunctionCalculator(leads(idx), 'shiftLead', shiftLeads(leads(idx)), 'transversepotential', transversepotential, 'coordinates_shift', coordinates_shift(leads(idx)),...
            'gauge_field', gauge_field, 'createCore', createCore, ...
            'SelfEnergy', SelfEnergy, 'SurfaceGreensFunction', SurfaceGreensFunction, 'q', q, 'leadmodel', leadmodel, 'CustomHamiltonian', CustomHamiltonian, 'Just_Create_Hamiltonians', Just_Create_Hamiltonians);
        obj.display(' ')
        idx = idx + 1;
    end

    lead_ret = obj.Leads;

end

%% DysonEq
%> @brief Invokes the function handle of the Dyson equation, and stores the calculated values in attributes #G and #junction_sites. 
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'CustomDyson' A function handle [G, Ginverz, #junction_sites] = CustomDyson() of the Dyson equation to be evaluated. If not given, the function handle given by the present class is used instead.
%> @return [1] Returns with the calculated Green operator.
%> @return [2] With the attribute #junction_sites.
function [Gret, junction_sites_ret] = DysonEq( obj, varargin )

    obj.display(['EQuUs:', class(obj), 'DysonEq: Evaluating Dyson Equation'])

    p = inputParser;
    p.addParameter('CustomDyson', []);        
    p.parse(varargin{:});

    if ~isempty( p.Results.CustomDyson )
        obj.Dysonfunc = p.Results.CustomDyson;
    end
    
    if ~isempty( obj.Dysonfunc )
        if nargout( obj.Dysonfunc ) == 1
            obj.G = obj.Dysonfunc();
        else
            [obj.G, ~, obj.junction_sites] = obj.Dysonfunc();            
        end
    else
        obj.G = [];
    end
    Gret = obj.G;
    junction_sites_ret = obj.junction_sites;

end


%% Conductance2
%> @brief Conductance calculated by Eq (19) in PRB 73 085414
%> @return Returns with the calculated conductance.
    function C = Conductance2( obj )
        
               
        Gamma = cell(size(obj.Leads));
        for idx = 1:length(obj.Leads)
            Gamma{idx} = obj.Leads{idx}.Gamma();
        end  
        
        Gamma_left = [Gamma{1}, zeros( size(Gamma{1},1), size(Gamma{2},2));
                      zeros( size(Gamma{2},1), size(Gamma{1},2)+size(Gamma{2},2)) ];
                  
        Gamma_right = [zeros( size(Gamma{1},1), size(Gamma{1},2)+size(Gamma{2},2));
                       zeros(size(Gamma{2},1), size(Gamma{1},2)), Gamma{2}]; 
                   
        for idx = 1:length(obj.Leads)
            Gamma{idx} = [];
        end
        
        % Eq (19) PRB 73 085414
        C =  trace( Gamma_left*obj.G'*Gamma_right*obj.G ) ;
        
        
%***********************************
%         Sigma = cell(size(obj.Leads));
%         for idx = 1:length(obj.Leads)
%             Sigma{idx} = obj.Leads{idx}.Read('Sigma');
%         end  
%         
%         Sigma_left = [Sigma{1}, zeros( size(Sigma{1},1), size(Sigma{2},2));
%                       zeros( size(Sigma{2},1), size(Sigma{1},2)+size(Sigma{2},2)) ];
%                   
%         Sigma_right = [zeros( size(Sigma{1},1), size(Sigma{1},2)+size(Sigma{2},2));
%                        zeros(size(Sigma{2},1), size(Sigma{1},2)), Sigma{2}]; 
%                    
%         for idx = 1:length(obj.Leads)
%             Sigma{idx} = [];
%         end
%         
%         Eq (22) Eur. Phys. J. B 53, 537-549 (2006)
%         current2 = 2*real( trace(-1i*obj.G*Gamma_right*obj.G'*Sigma_left') );
%         
%         current-current2
%*********************************
        

    end

%% Conduktance
%> @brief Calculates the conductance matrix from the scattering matrix
%> @return Returns with the calculated conductance matrix.
function C = Conduktance( obj )
    obj.display('calculating conductance matrix')

    obj.nyitott_csatornak = obj.nyitott_csatornak;
    LeadsNumber        = length(obj.open_channels);
    C = zeros(LeadsNumber, LeadsNumber);
    
    rows = 0;
    for rowidx = 1:LeadsNumber
        rows_tmp = max(rows)+ (1:obj.open_channels(rowidx).num);
        if isempty(rows_tmp)
            continue;
        end
        rows = rows_tmp;
        cols = 0;
        for colidx = 1:LeadsNumber
            cols_tmp = max(cols)+ (1:obj.open_channels(colidx).num);
            if isempty(cols_tmp)
                continue;
            end 
            cols = cols_tmp;
            s = obj.Smatrix( rows, cols);
           
            
            if ~obj.Opt.BdG
                if rowidx == colidx
                    C(rowidx,colidx) = -obj.open_channels(rowidx).num + trace(s'*s);
                else
                    C(rowidx,colidx) = trace(s'*s);      
                end
            else
                if obj.open_channels(rowidx).isSuperconducting || obj.open_channels(colidx).isSuperconducting
                    continue
                end                
                incoming_electrons = obj.open_channels(colidx).BdG_u_incoming( obj.open_channels(colidx).open_channels_p );
                outgoing_holes = ~obj.open_channels(rowidx).BdG_u_outgoing( obj.open_channels(colidx).open_channels_p );
                outgoing_electrons = obj.open_channels(rowidx).BdG_u_outgoing( obj.open_channels(colidx).open_channels_m );

                if rowidx == colidx
                    r_normal = s( outgoing_electrons,  incoming_electrons );
                    r_andreev = s( outgoing_holes,  incoming_electrons );
                    M = length(find( obj.open_channels(rowidx).open_channels_p & obj.open_channels(rowidx).BdG_u_incoming));
                    C(rowidx,colidx) = -M + trace(r_normal'*r_normal) - trace(r_andreev'*r_andreev);
                else
                    t_normal = s( outgoing_electrons,  incoming_electrons );
                    t_andreev = s( outgoing_holes,  incoming_electrons );                  
                    C(rowidx,colidx) =  trace(t_normal'*t_normal) - trace(t_andreev'*t_andreev);
                end                

                
            end
        end
    end
    
    obj.Conductance_Matrix = C;


end

%% SmatrixCal
%> @brief Calculates the scattering matrix
%> @return [1] The scattering matrix
%> @return [2] An array containing the number of openchannels in the leads. Obsolete, use attribute #open_channels istead.
function [S, Neff] = SmatrixCalc( obj )
    obj.display('calculating S-matrix')
    Neff = obj.Get_Neff();
    obj.Determine_Open_Channels()        

    NormamtxAll    = zeros(sum(Neff), sum(Neff));
    Modusmtx_pAll  = zeros(sum(Neff), sum(Neff));
    d_Modusmtx_mAll = zeros(sum(Neff), sum(Neff));
    Gyokvcsop_pAll = zeros(sum(Neff), 1);
    Gyokvcsop_mAll = zeros(sum(Neff), 1);
    evanescent_indexes_pAll = false(sum(Neff), 1);
    evanescent_indexes_mAll = false(sum(Neff), 1);

    for idx = 1:length(Neff)
        NormamtxAll( sum(Neff(1:idx-1))+1:sum(Neff(1:idx)), sum(Neff(1:idx-1))+1:sum(Neff(1:idx)) ) = obj.Leads{idx}.Read('Normamtx');
        Modusmtx_pAll( sum(Neff(1:idx-1))+1:sum(Neff(1:idx)), sum(Neff(1:idx-1))+1:sum(Neff(1:idx)) ) = obj.Leads{idx}.Read('modusmtx_p');
        d_Modusmtx_mAll( sum(Neff(1:idx-1))+1:sum(Neff(1:idx)), sum(Neff(1:idx-1))+1:sum(Neff(1:idx)) ) = obj.Leads{idx}.Read('d_modusmtx_m');
        Gyokvcsop_pAll( sum(Neff(1:idx-1))+1:sum(Neff(1:idx)) ) = sqrt(abs(obj.Leads{idx}.Read('vcsop_p')));
        Gyokvcsop_mAll( sum(Neff(1:idx-1))+1:sum(Neff(1:idx)) ) = sqrt(abs(obj.Leads{idx}.Read('vcsop_m')));
        evanescent_indexes_pAll( sum(Neff(1:idx-1))+1:sum(Neff(1:idx)) ) = ~(obj.open_channels(idx).open_channels_p);
        evanescent_indexes_mAll( sum(Neff(1:idx-1))+1:sum(Neff(1:idx)) ) = ~(obj.open_channels(idx).open_channels_m);
    end
    
    GN0 = obj.G(1:sum(Neff),1:sum(Neff))*NormamtxAll;
    S   = (d_Modusmtx_mAll*(GN0-eye(sum(Neff))) )*Modusmtx_pAll;

    % sorting out evanescent modes
    S(:,evanescent_indexes_pAll) = [];
    S(evanescent_indexes_mAll,:) = []; 
    Gyokvcsop_mAll = Gyokvcsop_mAll(~evanescent_indexes_mAll);
    Gyokvcsop_pAll = Gyokvcsop_pAll(~evanescent_indexes_pAll);
    
    S = diag( Gyokvcsop_mAll )*S*diag(Gyokvcsop_pAll.^(-1));  
    obj.Smatrix = S;

    % obsolete code lines
    obj.nyitott_csatornak = zeros(size(Neff));
    for idx = 1:length(Neff)      
        obj.nyitott_csatornak(idx) = obj.open_channels(idx).num;
    end
    Neff = obj.nyitott_csatornak;
 

end


%% Determine_Open_Channels
%> @brief Determines the open channels in the leads. The data are storen within the attribute #open_channels.
function Determine_Open_Channels( obj )
    M = obj.GetM();
    obj.open_channels = structures('open_channels');
    obj.open_channels = repmat(obj.open_channels, size(M));
    
    for idx = 1:length( M )
        obj.Leads{idx}.Determine_Open_Channels();
        obj.open_channels(idx) = obj.Leads{idx}.Read('open_channels');
    end
    
end

%% SurfaceGreenFunctionCalculator
%> @brief Calculates the surface Green's function or the self energy of a Lead.
%> @param idx the identification number of the lead in the attribute Leads (see attribute #CreateLeadHamiltonians.Hanyadik_Lead).
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'shiftLeads' A real number. If given, the on-site potentials of the sites in the lead are shifted by this value.
%> @param 'transversepotential' A function handle pot = f( #coordinates ) or pot=f( #CreateLeadHamiltonians, Energy) of the transverse potential applied in the lead. (Instead of #CreateLeadHamiltonians can be used its derived class)
%> @param 'Lead' An instance of class #Lead (or its subclass) that is intended to be used in the calculations. 
%> @param 'coordinates_shift' An integer. If given, the coordinates of the sites in the lead are shifted by coordinates_shift*lattice vector.
%> @param 'gauge_field' Function handle S = f( x,y) of the gauge transformation on the Hamiltonians of the leads. (S is a N x 1 vector, where N is the number of the points given by the x and y coordinates.)
%> @param 'createCore' Set 1 for creating the class instances #Lead but omitting any further calculations, or false (default) to continue with the calculations.
%> @param 'SelfEnergy' Logical value. Set true to calculate the self energies of the leads. (default is false)
%> @param 'SurfaceGreensFunction' Logical value. Set true to calculate the surface Green functions of the leads. (default is true)
%> @param 'leads' Array of lead identification numbers for which the calculations should be performed. (Default is a list of all leads)
%> @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 'CustomHamiltonian' An instance of class #Custom_Hamiltonians or its derived class to create custom Hamiltonians in the leads.
%> @param 'Just_Create_Hamiltonians' Logical value. Set true to create the Hamiltonians of the leads, but stop further calculations. Default value is false.
%> @param 'q' The tranverse momentum for transverse computations.
%> @return Returns with the created #Lead instance.
function Lead_ret = SurfaceGreenFunctionCalculator( obj, idx, varargin)
    
    p = inputParser;
    p.addParameter('createCore', 0);
    p.addParameter('Just_Create_Hamiltonians', 0);
    p.addParameter('shiftLead', 0);
    p.addParameter('coordinates_shift', 0);
    p.addParameter('transversepotential',[]);
    p.addParameter('Lead',[]);
    p.addParameter('gauge_field', [] );% gauge field for performing gauge transformation
    p.addParameter('SelfEnergy',false); % set true to calculate the self energy of the semi-infinite lead
    p.addParameter('SurfaceGreensFunction', true );% set true to calculate the surface Greens function of the semi-infinite lead
    p.addParameter('leadmodel', []); %function handle for an individual physical model for the contacts
    p.addParameter('CustomHamiltonian', []); 
    p.addParameter('q', [] ) %transverse momentum
    p.parse(varargin{:});
    createCore            = p.Results.createCore;
    Just_Create_Hamiltonians = p.Results.Just_Create_Hamiltonians;
    shiftLead             = p.Results.shiftLead;
    coordinates_shift     = p.Results.coordinates_shift;
    transversepotential             = p.Results.transversepotential;
    Lead_ret           = p.Results.Lead;
    gauge_field           = p.Results.gauge_field; % gauge field for performing gauge transformation
    SelfEnergy            = p.Results.SelfEnergy;
    SurfaceGreensFunction = p.Results.SurfaceGreensFunction;
    q                     = p.Results.q;
    CustomHamiltonian     = p.Results.CustomHamiltonian;
    leadmodel             = p.Results.leadmodel;
    
    
    if ~isempty( Lead_ret )
        % do not create a new class instance if Lead_ret was given
        Lead_ret.Reset();                                    
                                        
    elseif ~isempty(idx) && ~isempty(obj.param.Leads{idx}.Lead_Orientation)     
        Lead_ret = Lead(obj.Opt, obj.param,...
                                            'Hanyadik_Lead', idx,....
                                            'Lead_Orientation', obj.param.Leads{idx}.Lead_Orientation, ...
                                            'q', q);
    else
        Lead_ret = Lead(obj.Opt, obj.param,...
                                            'Hanyadik_Lead', idx, ...
                                            'q', q);
    end
    
    
    if createCore
        return
    end
                                        
    if obj.Opt.Simple_Green_Function
        Lead_ret.SurfaceGreen_simple(obj.E);       
    elseif ~isempty(leadmodel)
        Lead_ret = leadmodel( idx, obj.E, 'createCore', createCore, ...
                            'Just_Create_Hamiltonians', Just_Create_Hamiltonians, ...
                            'shiftLead', shiftLead, ...
                            'coordinates_shift', coordinates_shift, ...
                            'transversepotential', transversepotential, ...
                            'Lead', Lead_ret, ...
                            'gauge_field', gauge_field, ...
                            'SelfEnergy', SelfEnergy, ...
                            'SurfaceGreensFunction', SurfaceGreensFunction, ...
                            'q', q);
    else
    
        % creating Hailtonians
        obj.display(['EQuUs:', class(obj), ':SurfaceGreenFunctionCalculator: Creating Hamiltonians for lead'])
        Lead_ret.CreateHamiltonians( 'CustomHamiltonian', CustomHamiltonian );
        Lead_ret.ShiftCoordinates( coordinates_shift );
        if obj.Opt.magnetic_field == 1
            obj.display(['EQuUs:', class(obj), ':SurfaceGreenFunctionCalculator: Applying magnetic field in lead'])
            
            % In superconducting lead one must not include nonzero magnetic
            % field.
            % Hamiltonians in transverse computations must remain
            % traslational invariant. In transverse computations a gauge
            % tranformation is performed on each lead-sample interface to
            % cancel the vector potential
            if ~Lead_ret.isSuperconducting() %&& isempty( Lead_tmp.Read('q') )
                obj.PeierlsTransform.PeierlsTransformLeads(Lead_ret);
            elseif ~isempty( gauge_field ) %&& isempty( Lead_tmp.Read('q') )
               obj.PeierlsTransform.gaugeTransformationOnLead( Lead_ret, gauge_field );
            end
        end
        
        if abs(shiftLead) > 1e-6
            Lead_ret.ShiftLead( shiftLead );
        end

        if ~isempty(transversepotential) && isempty(q) %In transverse computations no transverse potential can be applied
            coordinates = Lead_ret.Read('coordinates');
               if nargin( transversepotential ) == 1
                    potential2apply = transversepotential( coordinates );
              elseif nargin( transversepotential ) == 2
             potential2apply = transversepotential( Lead_ret, obj.E ); 
         else
             error('EQuUs:Transport_Interface:SurfaceGreenFunctionCalculator', 'To many input arguments in function handle PotInScatter');
         end

               if isprop(coordinates, 'BdG_u')
                    if size( potential2apply, 1) == 1 || size( potential2apply, 2) == 1
                     potential2apply(~coordinates.BdG_u) = -potential2apply(~coordinates.BdG_u);
                    else
                         potential2apply(~coordinates.BdG_u, ~coordinates.BdG_u) = -conj(potential2apply(~coordinates.BdG_u, ~coordinates.BdG_u));
                    end
               end
               Lead_ret.AddPotential( potential2apply );
        end    
        
        if Just_Create_Hamiltonians
            return;
        end
                
        % Solve the eigen problem
        obj.display(['EQuUs:', class(obj), ':SurfaceGreenFunctionCalculator: Eigenvalues for lead'])
        Lead_ret.TrukkosSajatertekek(obj.E);
        
        % group velocities
        obj.display(['EQuUs:', class(obj), ':SurfaceGreenFunctionCalculator: Group velocities for lead'])
        Lead_ret.Group_Velocity();
               
        
        % retarded surface Green operator
        if SurfaceGreensFunction
            Lead_ret.SurfaceGreenFunction();
        end
        
        % retarded SelfEnergy
        if SelfEnergy
            Lead_ret.SelfEnergy();
        end
        
    end
    
    
end

%% setEnergy
%> @brief Sets the energy to be used in the calculations and resets the calculated attributes.
%> @param Energy The energy value
    function setEnergy( obj, Energy )
        obj.E = Energy; 
       
        % reset the attributes in the CreateHamiltonians class
        if ~isempty( obj.CreateH )
            obj.CreateH.Reset();
        end
       
        % reset the lead attributes
        if ~isempty( obj.Leads )
            for idx = 1:length( obj.Leads )
                obj.Leads{idx}.Reset();
            end
        end
       
    end
    

%% Read
%> @brief Query for the value of an attribute in the class.
%> @param MemberName The name of the attribute.
%> @return Returns with the value of the attribute.
    function ret = Read(obj, MemberName)
        if strcmp(MemberName, 'M')
            ret = obj.GetM();    
            return
        else
           try
                ret = obj.(MemberName);
           catch
                error(['EQuUs:', class(obj), ':Read'], ['No property to read with name: ',  MemberName]);
            end           
        end
        
    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')
            obj.param = input;
            obj.Reset();
            return 
        elseif strcmp(MemberName, 'E')
            obj.setEnergy( input );        
        else
               try
               obj.(MemberName) = input;
               catch
                    error(['EQuUs:', class(obj), ':Read'], ['No property to write 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    
    
%% Reset
%> @brief Resets all elements in the class.
    function Reset(obj)
        
        if strcmpi( class(obj), 'Transport_Interface' )
            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') || strcmp(prop_name, 'E')
                    continue
                end
                obj.Clear( prop_name ); 
            end   
        end 
        
        obj.Initialize();      
          
    
    end
    

%% replaceLead
%> @brief Adds/replaces a lead to/in the system.
%> @param Lead_tmp An instance of class #Lead
%> @param idx Identification number of the lead. (see attribute #CreateLeadHamiltonians.Hanyadik_Lead).
    function replaceLead( obj, Lead_tmp, idx )
        
        if ~strcmpi( class(obj.junction), 'Lead' )        
            supClasses = superclasses(Lead_tmp);
            if sum( strcmp( supClasses, 'Lead' ) ) == 0
                error(['EQuUs:', class(obj), ':replaceLead'], 'Input Lead_tmp is not valid.');
            end
        end
        
        obj.Leads{ idx } = Lead_tmp;   
        obj.Opt.NofLeads = length(obj.Leads);
    end

%% LeadSave
%> @brief Saves the Hamiltonians of each lead
    function LeadSave( obj )        
       for idx = 1:length(obj.Leads)
          obj.Leads{idx}.SaveLeads();
       end 
    end

%% CreateClone
%> @brief Creates a clone of the present class.
%> @return Returns with the cloned object.
    function ret = CreateClone( obj )
        
        if isempty( obj.CreateH )
            CreateH_loc = [];
        else
            CreateH_loc = obj.CreateH.CreateClone();
        end
        
        if isempty( obj.PeierlsTransform )
            PeierlsTransform_loc = [];
        else
            PeierlsTransform_loc = obj.PeierlsTransform.CreateClone();
        end
        
        ret = Transport_Interface(obj.E, obj.Opt, obj.param, ...
            'CreateH', CreateH_loc, ...
            'PeierlsTransform', PeierlsTransform_loc );
        
        for idx = 1:length(obj.Leads)
            ret.replaceLead( obj.Leads{idx}.CreateClone(), idx );
        end      
        
    end   

%% GetM
%> @brief Gets the width of the leads (see attribute #CreateLeadHamiltonians.M).
%> @return Returns with an array of the lead widths
    function M = GetM( obj )
        M = zeros(length(obj.Leads),1);
        for idx = 1:length(obj.Leads)
           M(idx) = obj.Leads{idx}.Read( 'M' );            
        end  
        
%         if ~isempty(M)
%             param.Leads.M = M;
%         end
    end
    
%% Get_Neff
%> @brief Gets the size of the effective Hamiltonians of the leads
%> @return Returns with an array of integers;
    function Neffs = Get_Neff( obj )
        Neffs = zeros(length(obj.Leads),1);
        for idx = 1:length(obj.Leads)
           Neffs(idx) = obj.Leads{idx}.Get_Neff();            
        end  
    end    
            

end % method public


methods ( Access = protected )
    
%% Initialize
%> @brief Initializes object attributes.
    function Initialize(obj)
                
        obj.G = [];
        obj.junction_sites = [];
        obj.Smatrix = [];
        obj.nyitott_csatornak  = [];
        obj.open_channels = [];
        obj.Conductance_Matrix = [];
        obj.Dysonfunc = [];

        obj.CreateH = [];
        obj.PeierlsTransform = [];
        obj.Leads = [];
        
        if strcmpi( class(obj), 'Transport_Interface')
            obj.InputParsing( obj.varargin{:});
        end
        
    end        

%% 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 'CreateH' An instance of class #CreateHamiltonians or its subclass storing the Hamiltonian of the scattering center.
%> @param 'PeierlsTransform' An instance of class #Peierls for the Peierls transformations.
    function InputParsing( obj, varargin )
        
        p = inputParser;
        p.addParameter('CreateH', []);
        p.addParameter('PeierlsTransform', []);
        
        p.parse(varargin{:});

        obj.PeierlsTransform = p.Results.PeierlsTransform;
        obj.CreateH = p.Results.CreateH;

    end

end % methods protected

end