antidot.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - antidot
%    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 utilities Utilities
%> @{
%> @file antidot.m
%> @brief A class to perform transport calculations on a graphene antidot (i.e., a hollow in a ribbon). Obsolete class, for real calculations a creation of a new class is recommended.
%> @} 
%> @brief A class to perform transport calculations on a graphene antidot (i.e., a hollow in a ribbon). Obsolete class, for real calculations a creation of a new class is recommended.
%%
classdef antidot < Ribbon

properties ( Access = protected )
    %> A lattice vector in the hexagonal lattice.
    a1
    %> A lattice vector in the hexagonal lattice.
    a2
    %> Directory name to export the plotted wave functions.
    waveFncDirnameFull
    %> The name of the subdirectory to export the plotted wave functions.
    waveFncSubDirname
end

properties ( Access = public )
    %> An instance of structure hole
    hole
    %> an instance of structures scatterers
    scatterers
    %> The strength of the magnetic field in Tesla
    B
    %> an instance of structure structures
    coordinates
    %> the sites at the edge of the antidot
    antidot_edge_points
    %> flux in the hole in units of \phi_0
    flux_of_hole
    %> atom-atom distance
    rCC
    %> Planck contant
    h = 6.626e-34;
    %> The charge of the electron
    qe = 1.602e-19;        
end



methods ( Access = public )
    
%% antidot
%> @brief Constructor of the class.
%> @param Opt An instance of the structure #Opt.
%> @param interpq An array of transverse momentum points at which fhandle is calculated via interpolation.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
    function obj = antidot( varargin )
        % Racsvektorok
        obj.a1 = [sqrt(3); -1]*sqrt(3)/2;
        obj.a2 = [sqrt(3); 1]*sqrt(3)/2;
        
        obj.hole = structures('hole');
        obj.scatterers = structures('scatterers');
        obj.coordinates = structures('coordinates');
            
        obj.B = 10e-10; %in Tesla
        obj.antidot_edge_points = [];                
        obj.flux_of_hole = [];
        

        obj.G      = [];
        obj.Ginv  = [];

        obj.waveFncDirnameFull = [];
        obj.waveFncSubDirname = [];

        obj.InputParsing( varargin{:} );
        obj.generate_geometry();
        
        obj.display('Creating antidot object')
        obj.createShape();
        obj.setFermiEnergy();
        
        createOutput( obj.filenameOut, obj.Opt, obj.param );
        
        obj.CreateHandles();
        % create Hamiltonians and coordinates
        obj.CreateRibbon('justHamiltonians', true);
    end

%% Transport
%> @brief Calculates the conductance of an antidot connected to the leads in the "contact/scattering center/contact" arrangement.
%> @param Energy The energy value.
%> @return C The calculated conductivity
%> @return ny Array of the open channel in the leads.
%> @return DeltaC Error of the unitarity.
    function [C,ny,DeltaC] = Transport( obj, Energy )
     obj.display(' ')
     obj.display( 'Calculating transport in create_Hole_in_ribbon')
        obj.setEnergy( Energy );
        obj.create_scatter_GreensFunction();
        
        Dysonfunc = @()obj.CustomDysonFunc('gfininv', obj.Ginv, 'SelfEnergy', false, 'constant_channels', false);
        try
            obj.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
        catch errCause
            save('Error_antidot_Transport.mat');
            for idx = 1:length(errCause.stack)
                obj.display(errCause.stack(idx), 1)
            end
            throw( errCause )
        end
        [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
        
        if ny(1) ~= ny(2)
           obj.display( ['openchannels do not match: ',num2str(ny)] )
        end
                
        conductance = obj.FL_handles.Conduktance();
        conductance = abs([conductance(1,:), conductance(1,:)]);
        DeltaC = std(conductance);
        
        C = mean(conductance);
        
        obj.display( ['Energy = ', num2str(Energy), ' conductance = ', num2str(C)])
        
        
    end

%% CalcWavefnc
%> @brief Calculates the energies and wave functions of the bound states.
%> @param Energy The energy value.
%> @param closefigure If true, the created figure is closed at the end of the calculations. 
%> @param varargin Cell array of optional parameters:
%> @param 'toPlot' Set true (default) for plotting the wave function, or false otherwise.
%> @param 'Einhomegac' Set true if the Energy is given in units of $$\hbar\omega_c$$, or false (default) otherwise.
%> @param 'filterHole' Set true (default) for separating antidot bound states from the edge states, or false otherwise.
%> @param 'db' The number of the calculated energy eigenvalues. Default is db=30.
%> @param 'infinitemass' Set true for using the infinite mass boundary condition, or false (default) otherwise.
%> @param 'dotCalc' Set true for a dot, or false (default) otherwise.
%> @param 'delta' Every delta-th point of the wave function becomes plotted. Default is delta=4.
%> @param 'smoothedge' Set true for using a smooth edge in the calculations, or false (default) otherwise.
%> @return ret A structure containing the calculated results.
%TODO create structure corresponding to the return value
    function ret = CalcWavefnc( obj, Energy, closefigure, varargin )
     obj.display( 'Calculating the wave function' );
        if ~exist('closefigure', 'var')
            closefigure = 1;
        end
    
        p = inputParser;
        p.addParameter('toPlot', 1, @isscalar);  
        p.addParameter('Einhomegac', 0, @isscalar);
        p.addParameter('filterHole', 1, @isscalar);
        p.addParameter('db', 30, @isscalar);
        p.addParameter('infinitemass', 0, @isscalar);
        p.addParameter('dotCalc', 0, @isscalar);
        p.addParameter('delta', 4, @isscalar);
        p.addParameter('smoothedge', 0, @isscalar);
        
        p.parse(varargin{:});    
        toPlot     = p.Results.toPlot;
        Einhomegac = p.Results.Einhomegac;
        filterHole = p.Results.filterHole; 
        db         = p.Results.db;
        infinitemass = p.Results.infinitemass;
        dotCalc      = p.Results.dotCalc;
        delta        = p.Results.delta; %used in wavefunction plotting
        smoothedge   = p.Results.smoothedge;
        
        if Einhomegac
            rCC = obj.rCC*1e-10; %Angstrom
            phi0 = obj.h/obj.qe; % magnetic flux quantum
            eta_B = 2*pi/phi0*(rCC)^2*obj.B;
            homega = 2.97*sqrt(9/2*eta_B);
            Energy = Energy*homega;
        end        
        
        
        
         if dotCalc && ~infinitemass
            filterHole = 0;
         end
            [Hscatter, wavefunc_indexes] = obj.create_hole_Hamiltonian( 'infinitemass', infinitemass, 'dotCalc', dotCalc, 'smoothedge', smoothedge );
            xcoords = obj.coordinates.x';
            ycoords = obj.coordinates.y';
        
            M = obj.width;
            height = obj.height*2;
                  
            obj.display(' Calculating eigenvalues and wave functions')
            [eigvecs,eigvals] = eigs(Hscatter,db,Energy);
            
            eigvals = diag( eigvals );
            [eigvals, IX] = sort(eigvals);
            eigvecs = eigvecs(:,IX);
            ret = struct( 'eigenval',[], 'Wavefnc', [], 'xcoords', [], 'ycoords', []);
            
            if filterHole
                filterHoleStates()
                ret.hole = 1;
            end
            obj.display( eigvals)
            ret.eigenval = eigvals;
            
            % reshaping the eigenvectors
            Wavefnc = cell( size(eigvals ) );
            for kdx = 1:length(eigvals)
         
                eigvec = eigvecs(:,kdx);         
                wavefnc = zeros(size(xcoords));              
            
                % plotting wavefunction calculated by the faster way
                sorok = round( (wavefunc_indexes - mod(wavefunc_indexes,M))/M ) + 1;
                indexes_enrows = logical( mod(wavefunc_indexes,M) == 0 );
                sorok( indexes_enrows ) = sorok( indexes_enrows ) -1;
                indexek_tmp = 0;
                for idx = 1:height
                    indexes_in_row = logical( sorok == idx );
                    indexes_in_row = wavefunc_indexes(indexes_in_row) - M*(idx-1);
                    wavefnc(idx, indexes_in_row ) = eigvec( indexek_tmp + (1:length(indexes_in_row) ) );
                    indexek_tmp = indexek_tmp + length(indexes_in_row);         
                end 
                
                Wavefnc{kdx} = wavefnc;
            end
            
            ret.eigenval = eigvals;
            ret.Wavefnc = Wavefnc;
            ret.xcoords = xcoords;
            ret.ycoords = ycoords;
        
         
         
         if ~toPlot
             return
         end
         
         mkdir( obj.waveFncDirnameFull );
         mkdir( obj.waveFncDirnameFull, obj.waveFncSubDirname);
         obj.display( ' Plotting wavefunctions' )
         for kdx = 1:length(eigvals)
            
            
                figure1 = figure('XVisual',...
    '0x21 (TrueColor, depth 24, RGB mask 0xff0000 0xff00 0x00ff)',...
    'Renderer','painters',...
    'Colormap',[1 1 1;0.948412716388702 0.948412716388702 0.948412716388702;0.896825432777405 0.896825432777405 0.896825432777405;0.845238089561462 0.845238089561462 0.845238089561462;0.793650805950165 0.793650805950165 0.793650805950165;0.742063522338867 0.742063522338867 0.742063522338867;0.690476179122925 0.690476179122925 0.690476179122925;0.638888895511627 0.638888895511627 0.638888895511627;0.58730161190033 0.58730161190033 0.58730161190033;0.581477105617523 0.581477105617523 0.581477105617523;0.575652658939362 0.575652658939362 0.575652658939362;0.569828152656555 0.569828152656555 0.569828152656555;0.564003705978394 0.564003705978394 0.564003705978394;0.558179199695587 0.558179199695587 0.558179199695587;0.552354753017426 0.552354753017426 0.552354753017426;0.546530246734619 0.546530246734619 0.546530246734619;0.540705800056458 0.540705800056458 0.540705800056458;0.534881293773651 0.534881293773651 0.534881293773651;0.52905684709549 0.52905684709549 0.52905684709549;0.523232340812683 0.523232340812683 0.523232340812683;0.517407834529877 0.517407834529877 0.517407834529877;0.511583387851715 0.511583387851715 0.511583387851715;0.505758881568909 0.505758881568909 0.505758881568909;0.499934434890747 0.499934434890747 0.499934434890747;0.494109958410263 0.494109958410263 0.494109958410263;0.488285452127457 0.488285452127457 0.488285452127457;0.482460975646973 0.482460975646973 0.482460975646973;0.476636499166489 0.476636499166489 0.476636499166489;0.470812022686005 0.470812022686005 0.470812022686005;0.464987546205521 0.464987546205521 0.464987546205521;0.459163069725037 0.459163069725037 0.459163069725037;0.445249050855637 0.445249050855637 0.445249050855637;0.431335002183914 0.431335002183914 0.431335002183914;0.417420983314514 0.417420983314514 0.417420983314514;0.403506934642792 0.403506934642792 0.403506934642792;0.389592915773392 0.389592915773392 0.389592915773392;0.375678867101669 0.375678867101669 0.375678867101669;0.361764848232269 0.361764848232269 0.361764848232269;0.347850799560547 0.347850799560547 0.347850799560547;0.333936780691147 0.333936780691147 0.333936780691147;0.320022732019424 0.320022732019424 0.320022732019424;0.306108713150024 0.306108713150024 0.306108713150024;0.292194694280624 0.292194694280624 0.292194694280624;0.278280645608902 0.278280645608902 0.278280645608902;0.264366626739502 0.264366626739502 0.264366626739502;0.25045257806778 0.25045257806778 0.25045257806778;0.236538544297218 0.236538544297218 0.236538544297218;0.222624525427818 0.222624525427818 0.222624525427818;0.208710491657257 0.208710491657257 0.208710491657257;0.194796457886696 0.194796457886696 0.194796457886696;0.180882424116135 0.180882424116135 0.180882424116135;0.166968390345573 0.166968390345573 0.166968390345573;0.153054356575012 0.153054356575012 0.153054356575012;0.139140322804451 0.139140322804451 0.139140322804451;0.12522628903389 0.12522628903389 0.12522628903389;0.111312262713909 0.111312262713909 0.111312262713909;0.0973982289433479 0.0973982289433479 0.0973982289433479;0.0834841951727867 0.0834841951727867 0.0834841951727867;0.0695701614022255 0.0695701614022255 0.0695701614022255;0.0556561313569546 0.0556561313569546 0.0556561313569546;0.0417420975863934 0.0417420975863934 0.0417420975863934;0.0278280656784773 0.0278280656784773 0.0278280656784773;0.0139140328392386 0.0139140328392386 0.0139140328392386;0 0 0]);

                

            % Create axes
            axes1 = axes('Parent',figure1, 'CLim',[0 0.028]);
            hold on
              
            contourf( xcoords(1:delta:end,1:delta:end), ycoords(1:delta:end,1:delta:end), abs(wavefnc(1:delta:end,1:delta:end)), 50, 'LineStyle', 'none','Parent', axes1 ); 
            axis equal
        
            text(30,30, ['E = ', num2str(eigvals(kdx)*1000),' meV']);
            set(gcf,'Renderer','painters')
                
            filename = ['wavefunc_width',num2str(obj.width),'_height',num2str(obj.height),'_radius',num2str(obj.hole.radius),'_B',num2str(obj.B),'_E',num2str(eigvals(kdx),'%6.5f'),'.eps'];
            if infinitemass
                filename =  ['infmass_',filename];
            end
            if smoothedge
                filename =  ['smoothedge_',filename];
            end
            if dotCalc
                filename = ['dot_', filename];
            end
            print('-depsc2', [obj.waveFncDirnameFull, filesep, obj.waveFncSubDirname,'/',filename]);
            

            if closefigure
             close(figure1);
            end
            
            
         
         end
         
         
         %----------------------------------------------
         % nested functions         
         
         %----------------------------------
         %this filters out eigenstates related to bound states to the Hole
         function filterHoleStates()
%             disp( ' Filtering out hole states' )
              hole_indexes = true(1,length(eigvals));
              
              tol = 1e-2;
              for iidx = 1:length(hole_indexes)
                 if norm(eigvecs(1:M-20,iidx)) > tol
                     hole_indexes(iidx) = false;
                 end                
              end
              
              eigvals = eigvals( hole_indexes );
              eigvecs = eigvecs(:, hole_indexes );
         end
         
         % end nested functions
         
         
    end

%% create_hole_Hamiltonian
%> @brief Creates the Hamiltonian for the antidot/dot.
%> @param varargin Cell array of optional parameters:
%> @param 'infinitemass' Set true for using the infinite mass boundary condition, or false (default) otherwise.
%> @param 'dotCalc' Set true for a dot, or false (default) otherwise.
%> @param 'smoothedge' Set true for using a smooth edge in the calculations, or false (default) otherwise.
%> @return Hscatter The matrix representation of the created Hamiltonian. 
%> @return wavefunction_indexes The list of the sites included in the antidot/dot system. 
    function [Hscatter, wavefunction_indexes] = create_hole_Hamiltonian( obj, varargin )
        p = inputParser;
        p.addParameter('infinitemass', 0, @isscalar);
        p.addParameter('dotCalc', 0, @isscalar);
        p.addParameter('smoothedge', 0, @isscalar);
        p.parse(varargin{:}); 
        
        infinitemass = p.Results.infinitemass;
        dotCalc      = p.Results.dotCalc;
        smoothedge   = p.Results.smoothedge;        
        
        obj.CreateHandlesForMagneticField(); 
        obj.setEnergy([]);
        
        obj.display(' Creating the Hamiltonian of a finite ribbon')
        obj.CreateScatter( );
        Hscatter = obj.CreateH.Read( 'Hscatter' );
        
        obj.display(' Creating hole in the ribbon ')
        
        
        if infinitemass
           Hscatter_tmp = Hscatter;
           V = -2;
        elseif dotCalc
            Hscatter_tmp = Hscatter;
            V = [];
        else
            Hscatter_tmp = Hscatter;
            V = [];
        end
        
        
        if isempty(obj.antidot_edge_points)
            obj.getHolePoints();
        end 
            
          obj.display( ' cutting out hole points ')
        M = obj.width;
        indexes2clear = [];
        for idx = size(obj.antidot_edge_points,1):-1:1
            indexek_tmp = (obj.antidot_edge_points(idx)-1)*(2*M) + (obj.antidot_edge_points(idx,2) : obj.antidot_edge_points(idx,3)); 
            Hscatter( indexek_tmp,: ) = 0;
         Hscatter( :, indexek_tmp ) = 0;
            indexes2clear = [indexes2clear, indexek_tmp];
        end 
        
        indexes2clear_tmp = getHoleSitesFromTheOddRows( indexes2clear );
        indexes2clear = sort( [ indexes2clear, indexes2clear_tmp]);
        
        
        if ~infinitemass 
            indexes_tmp = 1:size(Hscatter_tmp,1);
            indexes_tmp( indexes2clear) = [];
            
            if dotCalc
               wavefunction_indexes = indexes2clear;
               Hscatter = Hscatter_tmp;
               Hscatter( indexes_tmp,: ) = [];
               Hscatter( :, indexes_tmp ) = [];
            else
                if smoothedge 
                    wavefunction_indexes = indexes_tmp;
                    Hscatter( indexes2clear,: ) = [];
                    Hscatter( :, indexes2clear ) = [];
                else
                    wavefunction_indexes = 1:size(Hscatter,1);
                end
            end            
        else
            if dotCalc
               indexes_tmp = 1:size(Hscatter_tmp,1);
               indexes_tmp( indexes2clear ) = [];
               indexes2clear = indexes_tmp;
            end
            ApplyInfiniteMass( indexes2clear )  
            wavefunction_indexes = 1:size(Hscatter_tmp,1);
        end

        %----------------------------------------------
        % nested functions

        %-----------------------------------------
        function ApplyInfiniteMass( indexek)
            Hscatter = Hscatter_tmp;
            
            sorok = round((indexek - mod(indexek,M))/M)+1;
            indexes_endrows = logical( mod(indexek,M) == 0 );
            sorok( indexes_endrows ) = sorok( indexes_endrows ) -1;
            
            ps_sorok = logical( mod(sorok,2) == 0 );                
            pt_sorok = ~ps_sorok;
            
            ps_sorok = indexek( ps_sorok );
            pt_sorok = indexek( pt_sorok );
            
            idx_A = [ ps_sorok( logical( mod(ps_sorok,2) == 0)), pt_sorok( logical( mod(pt_sorok,2) == 1)) ];
            idx_B = [ ps_sorok( logical( mod(ps_sorok,2) == 1)), pt_sorok( logical( mod(pt_sorok,2) == 0)) ];
            
            Hscatter = Hscatter + sparse( idx_A, idx_A, V, size(Hscatter,1), size(Hscatter,2)) + sparse( idx_B, idx_B, -V, size(Hscatter,1), size(Hscatter,2));
            
            
        end
        
        
        % ----------------------------------------
        function indexes2clear = getHoleSitesFromTheOddRows( indexek_tmp )
            
            indexes2clear = [];

            for iidx = 1:length( indexek_tmp )
                related_indexes = find(Hscatter_tmp(indexek_tmp(iidx),:));
                sorok = round((related_indexes - mod(related_indexes,M))/M)+1;
                indexes_endrows = logical( mod(related_indexes,M) == 0 );
                sorok( indexes_endrows ) = sorok( indexes_endrows ) -1;
                
                ps_sorok = logical( mod(sorok,2) == 0);
                related_indexes = related_indexes( ps_sorok );
                
                for iiidx = 1:length(related_indexes)
                    idx_tmp = find(Hscatter(related_indexes(iiidx),:));
                    idx_tmp2 = find( idx_tmp ~= related_indexes(iiidx) );
                    if length( idx_tmp2 ) <= 1;
                        if isempty( find(indexes2clear == related_indexes(iiidx),1) )
                            indexes2clear = [indexes2clear, related_indexes(iiidx)];
                        end
                    end
                end
                
            end
            
            
        end
        
        % end nested functions
        
        
    end
   
%% CreateHandlesForMagneticField
%> @brief Creates function handles of the vector potentials and apply the magnetic filed in the ribbon Hamiltonians.
    function CreateHandlesForMagneticField( obj )        
        hLandauy = obj.createVectorPotential( obj.B );        
        obj.setHandlesForMagneticField('scatter', hLandauy, 'lead', hLandauy );          
    end
    
end


methods ( Access = protected )
 
    %% createVectorPotential
%> @brief creates handle for vector potential
    function hLandauy = createVectorPotential( obj, B )
    
     rCC = obj.rCC*1e-10; %Angstrom
          phi0 = obj.h/obj.qe;
          
        eta_B = 2*pi/phi0*(rCC)^2*B;
        Aconst = 0.1; % constant vektor potential in the systems: stabilizes the numerical computation        
                
        hLandauy = @(x,y)([zeros(size(x,1),1),eta_B*x]);
        
              
        
        
    end
    
%% generate_geometry
%> @brief creates geometry data of the hole in the ribbon
    function generate_geometry( obj )
        
        obj.display(' ')
        obj.display(' Generating geometry in antidot ')
        obj.display( ' ' )
        
        % creating geometry with the 
        obj.convert_site_indexes( )
    end
    
end % methods private 

methods ( Access = public )

%% getHolePoints
%> @brief Determines the sites that should be cut off from the ribbon in order to create the hole.
%> @returns antidot_edge_points A matrix with culomns z (slab index), zpoints_min, zpoints_max (lower and upper bounds of the sites in the slab z).
    function antidot_edge_points = getHolePoints( obj )
        
        if isempty(obj.Scatter_UC)
            obj.CreateRibbon('justHamiltonians', true);
        end
        
        M = obj.Scatter_UC.Read( 'M' );
           
        deltax_ps = 2;  % 2n-1 es 2n kozotti tavolsag rCC egysegekben
        deltax_pt = 1;    % 2n es 2n+1 kozotti tavolsag rCC egysegekben
        deltax    = 3; % 2n-1 es 2n+1 kozotti tavolsag rCC egysegekben
        deltay = sqrt(3);            
            
        hole.center = [obj.hole.center(1)*(deltax_ps+deltax_pt)/2; obj.hole.center(2)*deltay];
            
        rA2 = obj.a2-obj.a1;
        hole.radius = obj.hole.radius*(norm(rA2))/2;
            
        zmin = round((hole.center(2) - hole.radius)/deltay) - 1;
        zmax = round((hole.center(2) + hole.radius)/deltay) + 1;
            
        antidot_edge_points = zeros( zmax-zmin, 3);
        hole_edge_idx = 1;
            
        for ztmp = zmin:zmax
            minfound = 0;
            for idx = 1:M
                dist = distance_from_center(ztmp, idx );
                if ( dist<hole.radius && ~minfound )
                    idx_min = idx;
                    minfound = 1;
                end
                    
                if ( dist>hole.radius && minfound )
                    idx_max = idx-1;
                    antidot_edge_points(hole_edge_idx,:) = [ztmp, idx_min, idx_max];
                    hole_edge_idx = hole_edge_idx + 1;
                    break;
                end
                    
            end  
        end
            
        antidot_edge_points(hole_edge_idx:end,:) = [];
        obj.antidot_edge_points = antidot_edge_points;
        
            %----------------------------------------------
            % nested functions
            
            %----------------------------------------------
            function ret = distance_from_center(z, idx )
                
                maradek = mod(idx,2);
                
                if maradek==1
                    xcoord = (idx-1)/2*deltax;
                else
                    xcoord = (idx-2)/2*deltax + deltax_ps;
                end
                
                ycoord = z*deltay;
                ret = sqrt( (xcoord-hole.center(1))^2 + (ycoord - hole.center(2))^2  );
                
            end
            
            % end nested functions
            
            
            
        end

%% create_scatter_GreensFunction        
%> @brief Calculates the surface Greens function of the antidot. Sets the attributes #gfin and #gfininv to the calculated results.
    function create_scatter_GreensFunction( obj )
        obj.CreateRibbon('justHamiltonians', true);
        obj.Scatter_UC.ApplyOverlapMatrices(obj.E);
        K0 = obj.Scatter_UC.Read('K0');
        K1 = obj.Scatter_UC.Read('K1');
        K1adj = obj.Scatter_UC.Read('K1adj');
        obj.CreateRibbon();
        
        if isempty(obj.antidot_edge_points)
            obj.getHolePoints();
        end        
        
        %> the first two and last two slabs are added manualy at the end
        z1 = 1;
        z2 = obj.height-2;
        
        
        obj.G = GreensFuncAtRibbonEnds();
        antidot_edge_points = obj.antidot_edge_points;
        scatterer_points = getScattererPoints();
        
        if ~isempty( scatterer_points )
            ghole = GreensFuncAtHoleEdge();
            [ghopp1, ghopp2] = GreensFuncAtCouplings();
            obj.G = [obj.G,ghopp1; ghopp2,ghole];
        end
        
        if isnan(rcond(obj.G) ) || abs( rcond(obj.G) ) < 1e-15
            obj.Ginv= obj.inv_SVD( obj.G );
        else
            obj.Ginv= inv(obj.G);%gfin\eye(size(gfin));%inv(gfin);
        end       
        obj.G = [];
        
        if ~obj.scatterers.aremissing
         %applying scattering potentials
            for ldx = 1:length(obj.scatterers.siteindexes)
                obj.Ginv( obj.scatterers.siteindexes(ldx)+4*M, obj.scatterers.siteindexes(ldx)+4*M) = ...
                 obj.Ginv( obj.scatterers.siteindexes(ldx)+4*M, obj.scatterers.siteindexes(ldx)+4*M) - obj.scatterers.potentials(ldx);
            end
            obj.Ginv= obj.Ginv( [1+M:3*M,obj.scatterers.siteindexes+4*M], [1+M:3*M,obj.scatterers.siteindexes+4*M] );
            kulso_szabfokok = 1:2*M;
            obj.Ginv= obj.DecimationFunction( kulso_szabfokok, gfininv);
        else
            obj.Ginv= obj.Ginv( 1+M:3*M, 1+M:3*M );
        end 
        
        % adding to the scattering region the first and last slabs
        Neff = obj.Scatter_UC.Get_Neff();
        non_singular_sites = obj.Scatter_UC.Read( 'kulso_szabfokok' ); 
        if isempty( non_singular_sites )
            non_singular_sites = 1:Neff;            
        end
        
        % getting the Hamiltonian for the edge slabs
        non_singular_sites_edges = [non_singular_sites, size(K0,1)+1:2*size(K0,1)];        
        ginv_edges = -[K0, K1; K1adj, K0];
        ginv_edges = obj.DecimationFunction( non_singular_sites_edges, ginv_edges );
        
        %Ginv = [first_slab, H1, 0;
        %             H1',  invG, H1;
        %              0  , H1', last_slab];        
                        
        obj.Ginv = [ginv_edges(1:Neff,1:Neff), [ginv_edges(1:Neff, Neff+non_singular_sites), zeros(Neff, Neff+size(K0,2))]; ...
            [ginv_edges(Neff+non_singular_sites, 1:Neff); zeros(Neff, Neff)], obj.Ginv, [zeros(Neff, size(ginv_edges,2)-Neff); ginv_edges(1:Neff, Neff+1:end)]; ...
                            [zeros(size(K0,2),2*Neff), ginv_edges(Neff+1:end, 1:Neff)], ginv_edges(Neff+1:end, Neff+1:end)];                        
                        
        [rows, cols] = find( K1 );
        rows = unique(rows);   % cols identical to non_singular_sites

        %non_singular_sites_Ginv = [1:length(non_singular_sites), size(obj.Ginv,1)-size(K0,1)+1:size(obj.Ginv,1)];
        non_singular_sites_Ginv = [1:length(non_singular_sites), size(obj.Ginv,1)-size(K0,1)+reshape(rows, 1, length(rows))];        
        obj.Ginv = obj.DecimationFunction( non_singular_sites_Ginv, obj.Ginv );        
        
        % apply gauge transformation if given
        coordinates_scatter = obj.getCoordinates();
        if ~isempty( obj.gauge_field )
            try
               % gauge transformation on Green's function
                if ~isempty(obj.Ginv) && ~isempty(obj.PeierlsTransform_Scatter) && isempty(obj.q)
                    % gauge transformation on the inverse Green's function
                    obj.Ginv= obj.PeierlsTransform_Scatter.gaugeTransformation( obj.Ginv, coordinates_scatter, obj.gauge_field );
                end
            catch  errCause  
                err = MException('Ribbon:CalcFiniteGreensFunction', 'Unable to perform gauge transformation');
                err = addCause(err, errCause);
                save('Error_Ribbon_CalcFiniteGreensFunction.mat')
                throw(err);
            end
        end 
        
        %----------------------------------------------
        % nested functions
        
       
        %--------------------------------------------------------------------------
          % get sites of hole and scatterer that should be obtained in the Green's function
        function scatterer_points = getScattererPoints()
          
              scatterer_points = cell(size(antidot_edge_points,1),3);
              
            if ~isempty(antidot_edge_points)      
                   point_num_down = antidot_edge_points(1,3)-antidot_edge_points(1,2) + 1;
             point_num_up = antidot_edge_points(end,3)-antidot_edge_points(end,2) + 1;
            end
               
            for idx = 1:size( antidot_edge_points,1 )
                ztmp = antidot_edge_points( idx );                
                
                if idx == 1 
                    zpoints = antidot_edge_points(idx,2):antidot_edge_points(idx,3);
                    indexek = 1:point_num_down;
                elseif idx == size(antidot_edge_points,1)
                    zpoints = antidot_edge_points(idx,2):antidot_edge_points(idx,3);
                    indexek = max(indexek)+1:max(indexek)+point_num_up;
                else
                    if antidot_edge_points(idx,2) + 6 < antidot_edge_points(idx,3)
                        zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,2)+2,antidot_edge_points(idx,2)+3,antidot_edge_points(idx,3)-3,antidot_edge_points(idx,3)-2,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
                        indexek = max(indexek)+1:max(indexek)+8;
                    elseif antidot_edge_points(idx,2) + 5 < antidot_edge_points(idx,3)
                        zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,2)+2,antidot_edge_points(idx,2)+3,antidot_edge_points(idx,3)-2,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
                        indexek = max(indexek)+1:max(indexek)+7;
                    elseif antidot_edge_points(idx,2) + 4 < antidot_edge_points(idx,3)
                        zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,2)+2,antidot_edge_points(idx,3)-2,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
                        indexek = max(indexek)+1:max(indexek)+6;                        
                    elseif antidot_edge_points(idx,2) + 3 < antidot_edge_points(idx,3)
                        zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,2)+2,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
                        indexek = max(indexek)+1:max(indexek)+5;                        
                    elseif antidot_edge_points(idx,2) + 2 < antidot_edge_points(idx,3)
                        zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
                        indexek = max(indexek)+1:max(indexek)+4;
                    elseif antidot_edge_points(idx,2) + 1 < antidot_edge_points(idx,3)
                        zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,3)];
                        indexek = max(indexek)+1:max(indexek)+3;
                    else
                        zpoints = antidot_edge_points(idx,2:3);
                        indexek = max(indexek)+1:max(indexek)+2;
                    end                   
                                                           
                end  
                
                scatterer_points{idx,1} = antidot_edge_points(idx,1);
                scatterer_points{idx,2} = zpoints;    
                scatterer_points{idx,3} = ones(size(zpoints)); %1 for hole edge points, 0 otherwise
                
            end                  
            
                  
          z = antidot_edge_points(:,1);
            
            for idx = 1:length(  obj.scatterers.z )
                    index = find( z == round(obj.scatterers.z(idx)/2),1 );
               if isempty( index )      
                    is_hole_edge = 0; %1 for hole edge points, 0 otherwise
                         scatterer_points = [scatterer_points; {round(obj.scatterers.z(idx)/2), obj.scatterers.zpoints(idx), is_hole_edge} ];
                         z = [z; round(obj.scatterers.z(idx)/2)];
                    else
                         scatterer_points{index,2} = [scatterer_points{index,2}, obj.scatterers.zpoints(idx)];
                    is_hole_edge = 0; %1 for hole edge points, 0 otherwise
                    scatterer_points{index,3} = [scatterer_points{index,3}, is_hole_edge];
                end
            end   
            
            getScattererSiteIndexes();
            

            
            %-----------------------------------------------------------
            function getScattererSiteIndexes()
                obj.scatterers.siteindexes = zeros( 1,length(  obj.scatterers.z ) );    
                
                index_sum = 0;
                kkdx = 0;
                for iidx = 1:size(scatterer_points,1)
                    
                    scatterer_idx = find( scatterer_points{iidx,3}==0 );
                    if ~isempty( scatterer_idx )
                        obj.scatterers.siteindexes(kkdx+1:kkdx+length(scatterer_idx) ) = index_sum + scatterer_idx;
                        kkdx = kkdx + length(scatterer_idx);
                    end
                    index_sum = index_sum + length(scatterer_points{iidx,3});                    
                end  
            end
          
        end
            

        
        %-------------------------------------------------
        % calculate the Greens function for the hole edge points
        function ghole = GreensFuncAtHoleEdge( )
            
            point_num = 0;
            for idx = 1:size(scatterer_points,1)
               point_num = point_num + length(scatterer_points{idx,2});
            end
            
            ghole = zeros(point_num,point_num);
            indexek1 = 0;
            for idx = 1:size( scatterer_points,1 )
                z1tmp = scatterer_points{idx,1};
                z1points = scatterer_points{idx,2};
                indexek1 = max(indexek1)+1:max(indexek1)+length(z1points);
                
                indexek2 = 0;
                for jdx = 1:size( scatterer_points,1 )
                         z2tmp    = scatterer_points{jdx,1};
                    z2points = scatterer_points{jdx,2};
                    indexek2 = max(indexek2)+1:max(indexek2)+length(z2points);              
                    
                                        
                    gtmp = obj.Scatter_UC.InfGreenFunction(z1tmp, z2tmp, 'z1points', z1points, 'z2points', z2points);  
                    
                    ghole(indexek1, indexek2) = gtmp;
                end
            end
            
            ghole(max(indexek1)+1:end,:) = [];
            ghole(:,max(indexek2)+1:end) = [];
            
        end
        
        
        %-------------------------------------------------
        % calculate the Greens function between the ribbon and points and hole edge
        function [ghopp1, ghopp2] = GreensFuncAtCouplings()
                  
                       
            point_num = size( ghole,1 );
            M = obj.Scatter_UC.Read( 'M' );
            ghopp1 = zeros(4*M,point_num);
            ghopp2 = zeros(point_num,4*M);
            
            indexek = 0;
            for idx = 1:size( scatterer_points,1 )
                ztmp    = scatterer_points{idx,1};
                zpoints = scatterer_points{idx,2};
                indexek = max(indexek)+1:max(indexek)+length(zpoints);                                                    
                
                gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, ztmp, 'z2points', zpoints);                    
                ghopp1(1:M, indexek) = gtmp;
                
                gtmp = obj.Scatter_UC.InfGreenFunction(z1, ztmp, 'z2points', zpoints);                    
                ghopp1(1+M:2*M, indexek) = gtmp;
                
                gtmp = obj.Scatter_UC.InfGreenFunction(z2, ztmp, 'z2points', zpoints);                    
                ghopp1(1+2*M:3*M, indexek) = gtmp;
                
                gtmp = obj.Scatter_UC.InfGreenFunction(z2+1, ztmp, 'z2points', zpoints);                    
                ghopp1(1+3*M:4*M, indexek) = gtmp; 
                
                
                
                
                gtmp = obj.Scatter_UC.InfGreenFunction( ztmp, z1-1, 'z1points', zpoints);                    
                ghopp2(indexek,1:M) = gtmp;
                
                gtmp = obj.Scatter_UC.InfGreenFunction( ztmp, z1, 'z1points', zpoints);                    
                ghopp2(indexek,1+M:2*M) = gtmp;
                
                gtmp = obj.Scatter_UC.InfGreenFunction( ztmp, z2, 'z1points', zpoints);                    
                ghopp2(indexek,1+2*M:3*M) = gtmp;
                
                gtmp = obj.Scatter_UC.InfGreenFunction( ztmp, z2+1, 'z1points', zpoints);                    
                ghopp2(indexek,1+3*M:4*M) = gtmp;                
                
            end
        end        
                
        
        %-------------------------------------------------
        function gfin = GreensFuncAtRibbonEnds()
            
            M = obj.Scatter_UC.Read('M');
            gfin = zeros(4*M,4*M);
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, z1-1);
            gfin(1:M,1:M)              = gtmp;
            gfin(1+M:2*M, 1+M:2*M)     = gtmp;
            gfin(1+2*M:3*M, 1+2*M:3*M) = gtmp;              
            gfin(1+3*M:4*M, 1+3*M:4*M) = gtmp;
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, z1);
            gfin(1:M,1+M:2*M) = gtmp;
            gfin(1+2*M:3*M,1+3*M:4*M) = gtmp;
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z1, z1-1);
            gfin(1+M:2*M,1:M) = gtmp;
            gfin(1+3*M:4*M, 1+2*M:3*M) = gtmp;    
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, z2);
            gfin(1:M,1+2*M:3*M) = gtmp; 
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z2, z1-1);
            gfin(1+2*M:3*M,1:M) = gtmp; 
                
            gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, z2+1);
            gfin(1:M,1+3*M:4*M) = gtmp; 
                
            gtmp = obj.Scatter_UC.InfGreenFunction(z2+1, z1-1);
            gfin(1+3*M:4*M,1:M) = gtmp;   
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z1, z2);
            gfin(1+M:2*M,1+2*M:3*M) = gtmp; 
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z2, z1);
            gfin(1+2*M:3*M,1+M:2*M) = gtmp; 
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z1, z2+1);
            gfin(1+M:2*M,1+3*M:4*M) = gtmp; 
        
            gtmp = obj.Scatter_UC.InfGreenFunction(z2+1, z1);
            gfin(1+3*M:4*M,1+M:2*M) = gtmp;  
            
        end
        
        % end nested functions
        
        
        
    end


end % methods public

methods ( Access = private )
  
%% convert_site_indexes
%> @brief Converets site indexes into coordinates
    function convert_site_indexes( obj )
        
        xcoord = zeros( obj.width, obj.height*2 );
        ycoord = zeros( obj.width, obj.height*2 );
        
        
        
        rA1 = obj.a1+obj.a2;
        rA2 = obj.a2-obj.a1;
        
        radius = obj.hole.radius*(norm(rA2))/2;
        
        obj.flux_of_hole = pi*radius^2*(obj.rCC*1e-10)^2*obj.B/(obj.h/obj.qe);
        
        rB = [-1;0];
        
        
        for idx = 1:obj.width
        for jdx = 1:obj.height*2
            
            if mod(idx+jdx,2) == 0 % A atom
                r = double(idx-1)/2*rA1 + double(jdx-1)/2*rA2;
            else %B atom
                r = double(idx)/2*rA1 + double(jdx-1)/2*rA2 + rB;
            end
            
            xcoord(idx,jdx) = r(1);
            ycoord(idx,jdx) = r(2);           
               
               
        end
        end
        
        %center = obj.hole.center(1)/2*rA1 + obj.hole.center(2)/2*rA2;
        %obj.coordinates.outhole = logical( (xcoord - center(1)).^2 + (ycoord - center(2)).^2 >= radius^2 );        
        
        obj.coordinates.x = xcoord*obj.rCC;
        obj.coordinates.y = ycoord*obj.rCC;
             
        
        
        
    end
    
end % methods private
    
methods ( Access = protected )    

%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Cell array of optional parameters:
%> @param 'width' Integer. The number of the atomic sites in the cross section of the ribbon.
%> @param 'height' Integer. The height of the ribbon in units of the lattice vector.
%> @param 'filenameIn' Input filename for the xml input structure.
%> @param 'filenameOut' Output filename for the xml input structure.
%> @param 'WorkingDir' The absolute path to the working directoy.
%> @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. (overrides the one comming from the external source)
%> @param 'silent' Set true for suppress the output messages.
%> @param 'transversepotential' A function handle pot=f( #coordinates coords) to calculate the transverse potential in the cross section of the ribbon.
%> @param 'leadmodel' A function handle #Lead clead=f( idx, E, varargin ) of the alternative lead model with equivalent inputs and return values as class #Lead 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.)
%> @param 'Opt' An instance of the structure #Opt. (Overrides data in the input file if given)
%> @param 'param' An instance of the structure #param. (Overrides data in the input file if given)
%> @param 'q' The transverse momentum quantum number.
%> @param 'B' The strength of the magnetic field in the units of Tesla.
%> @param 'hole' An instance of structure #hole.
%> @param 'scatterers' An instance of structure #scatterers.
%> @param 'radius' The radius of the hole in units of lattice constant.
%> @param 'rCC' The atomic distance.
    function InputParsing(obj, varargin)
        
        p = inputParser;
        p.addParameter('width', obj.width);
        p.addParameter('height', obj.height);
        p.addParameter('filenameIn', [pwd, filesep, 'Basic_Input_zigzag_leads.xml']);
        p.addParameter('filenameOut', [pwd, filesep, 'Basic_Input_running_parameters.xml']);
        p.addParameter('WorkingDir', pwd);
        p.addParameter('E', obj.E);
        p.addParameter('EF', obj.EF); 
        p.addParameter('silent', obj.silent);   
        p.addParameter('transversepotential', obj.transversepotential);
        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', []);
        p.addParameter('param', []);
        p.addParameter('q', []);
        
        p.addParameter('B', obj.B);
        p.addParameter('hole', obj.hole);
        p.addParameter('scatterers', obj.scatterers);
        p.addParameter('radius', 0);  
        p.addParameter('rCC', 1.42);
        
        p.parse(varargin{:});

          InputParsing@Ribbon( obj, 'width', p.Results.width, ...
                                   'height', p.Results.height, ...
                            'filenameIn', p.Results.filenameIn, ...
                            'filenameOut', p.Results.filenameOut, ...
                            'WorkingDir', p.Results.WorkingDir, ...
                            'E', p.Results.E, ...
                            'EF', p.Results.EF, ...
                            'silent', p.Results.silent, ...
                            'transversepotential', p.Results.transversepotential, ...
                            'leadmodel', p.Results.leadmodel, ...
                            'interfacemodel', p.Results.interfacemodel, ...
                            'Opt', p.Results.Opt, ...
                            'param', p.Results.param, ...
                            'q', p.Results.q);
   
        
        obj.B               = p.Results.B;
        obj.hole          = p.Results.hole;
        obj.scatterers = p.Results.scatterers;  
        obj.rCC        = p.Results.rCC;                        

        if p.Results.radius > 0
            obj.hole.radius = p.Results.radius;
        end
        
        obj.waveFncDirnameFull = [pwd, filesep, 'wave_function_full'];
        obj.waveFncSubDirname = ['width',num2str(obj.width),'_height',num2str(obj.height),'_radius',num2str(obj.hole.radius),'_B',num2str(obj.B)];

        
        

    end

end % methos protected


end