CombineRibbons.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - CombineRibbons
%    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/.
%
%
%>   An object for combining multiple ribbon parts of equal width and from the same material in a two terminal setup.
%%
classdef CombineRibbons < Ribbon
    
properties ( Access = private )
    %> An instance of structure coordinates.
    coordinates
end


properties
    %>  the list of the ribbon interface to be combained (Ribbon objects)
    Ribbons = {};
    %> the list of the interface regions between the ribbons (Surface_Green_function objects)      
    Interface_Regions_all = {};          
end



methods ( Access = public )
%% Contructor of the class
%> @brief Constructor of the class.
%> @param Opt An instance of the structure Opt.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class    
    function obj = CombineRibbons( varargin )
        obj = obj@Ribbon();     
        
        obj.coordinates = structures('coordinates');

        obj.InputParsing( varargin{:} );  
    end


%% Transport
%> @brief Calculates the transport through the two terminal setup on two dimensional lattices. Use for development pupose only.
%> @param Energy The energy value.
%> @param varargin Cell array of optional parameters:.
%> @return Conductance The conductance tensor
%> @return ny Array of the open channel in the leads.
%> @return DeltaC Error of the unitarity.
%> @return S The scattering matrix.
    function [Conductance,ny,DeltaC] = Transport( obj, Energy )
        obj.setEnergy( Energy )
        
        obj.CalcFiniteGreensFunction('onlyGinv', true);
        
        Surface_sc = cell(2,1);
        ribbon_tmp = obj.Ribbons{1};
        Surface_sc{1} = ribbon_tmp.Scatter_UC.CreateClone();
        Surface_sc{1}.ShiftCoordinates( -1 )
        ribbon_tmp = obj.Ribbons{end};
        total_height = obj.getTotalHeight();
        Surface_sc{2} = ribbon_tmp.Scatter_UC.CreateClone();
        Surface_sc{2}.ShiftCoordinates( total_height )
        
        
        Dysonfunc = @()(obj.CustomDysonFunc( 'constant_channels', false, 'gfininv', obj.Ginv ));
        try
            obj.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
        catch errCause
            err = MException('CombineRibbons:Transport', 'Error occured in calculating the Greens function');
            err = addCause(err, errCause);
            save('Error_CombineRibbons_Transport.mat');
            Conductance = NaN;
            ny = NaN(2,1);
            DeltaC = NaN;
            throw(err);
        end
        [S,ny] = obj.FL_handles.SmatrixCalc();
        
        norma = norm(S*S'-eye(sum(ny)));
            
        if norma >= 1e-3
            disp( ['error of the unitarity of S-matrix: ',num2str(norma)] )
        end
        
        if ny(1) ~= ny(2)
           disp( ['openchannels do not match: ',num2str(ny)] )
        end
        
        
        conductance = obj.FL_handles.Conduktance();
        conductance = abs([conductance(1,:), conductance(2,:)]);
        DeltaC = std(conductance);
        
        C = mean(conductance);
        Conductance = C;
        
        disp( ['conductance = ', num2str(C), ' open_channels= ', num2str(ny(1)), ' DeltaC = ', num2str(DeltaC)])        
        %2/pi
        
    end

%% getCoordinates
%> @brief Gets the coordinates of the central region
%> @return coordinates Coordinates of the surface sites of the central region.
%> @return coordinates_interface Coordinates of the interface region.
function [coordinates, coordinates_interface] = getCoordinates( obj )
   if isempty( obj.Ribbons )
       coordinates = [];
       coordinates_interface = [];
       return
   end
   
   [coordinates, coordinates_interface] = obj.Ribbons{1}.getCoordinates();    
    
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 )
        
        setEnergy@Ribbon( obj, Energy );
        
        for idx = 1:length(obj.Ribbons)            
            ribbon_loc = obj.Ribbons{idx};
            ribbon_loc.setEnergy( Energy )
        end
        
        
    end



%% CalcFiniteGreensFunctionFromHamiltonian
%> @brief Calculates the Green operator of the scattering region from the whole Hamiltonian.
%> @param varargin Cell array of optional parameters:. 
    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('E', [], @isscalar); %The energy to be used in the ribbons
        p.addParameter('onlyGinv', false ); %
        p.addParameter('ContactPotInterface', []);
        p.addParameter('PotInScatter', []);
        p.parse(varargin{:});

        gauge_trans = p.Results.gauge_trans;
        E           = p.Results.E;
        onlyGinv    = p.Results.onlyGinv;
        ContactPotInterface  = p.Results.ContactPotInterface;
        PotInScatter = p.Results.PotInScatter;
    
        obj.CalcFiniteGreensFunction( 'E', E, 'gauge_trans', gauge_trans, 'onlyGinv', onlyGinv, 'ContactPotInterface', ContactPotInterface, 'gfininvfromHamiltonian', true, 'PotInScatter', PotInScatter )        
    end


%% CalcFiniteGreensFunction
%> @brief Calculates the Green operator of the scattering region by the fast way (see PRB 90, 125428 (2014)).
%> @param varargin Cell array of optional parameters:. 
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('E', []); %The energy to be used in the ribbons
    p.addParameter('onlyGinv', false ); 
    p.addParameter('ContactPotInterface', []);
    p.addParameter('gfininvfromHamiltonian', false); %Not added to the documentation
    p.addParameter('PotInScatter', []); 
     p.parse(varargin{:});

    gauge_trans = p.Results.gauge_trans;
    E           = p.Results.E;
    onlyGinv    = p.Results.onlyGinv;
    ContactPotInterface  = p.Results.ContactPotInterface;
    gfininvfromHamiltonian = p.Results.gfininvfromHamiltonian;
    PotInScatter = p.Results.PotInScatter;
    
    obj.Ginv = [];
     
    if ~isempty( E )
        obj.setEnergy( E );
        for idx = 1:length( obj.Ribbons )
            ribbon_tmp = obj.Ribbons{idx};        
               ribbon_tmp.setEnergy( E );
        end
    end
    
    % creating the interface regions if needed
    if ~isempty( obj.Interface_Regions )
       obj.createInterfaceRegions()        
    end
    
    for idx = 1:length(obj.Ribbons)
       obj.attachRibbon( obj.Ribbons{idx}, 'onlyGinv', onlyGinv, 'gfininvfromHamiltonian', gfininvfromHamiltonian, 'PotInScatter', PotInScatter );
    end
    
    
    
    % gauge transformation of the vector potential in the effective Hamiltonians
    if gauge_trans
        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('CombineRibbons:CalcFiniteGreensFunction', 'Unable to perform gauge transformation');
            err = addCause(err, errCause);
            save('Error_CombineRibbons_CalcFiniteGreensFunction.mat')
            throw(err);
        end
    end 
        
end


%% setHandlesForMagneticField
%> @brief Sets the function handles of the vector potential and gauge transformation.
%> @param varargin Cell array of optional parameters:. % 
    function setHandlesForMagneticField( obj, varargin )        
    
        setHandlesForMagneticField@Ribbon( obj, varargin{:} );
        
        % other Ribbon classes in the list
        for idx = 1:length(obj.Ribbons)
          obj.Ribbons{idx}.setHandlesForMagneticField('scatter', hscatter, 'lead', hscatter);%hlead, 'gauge_field', gauge_field );         
        end      
    end
    
    
%% CreateClone
%> @brief Creates a clone of the present object.
%> @return Returns with the cloned object.
    function ret = CreateClone( obj )
        Ribbons_cloned = cell(size(obj.Ribbons));
        for idx = 1:length(Ribbons_cloned)
           Ribbons_cloned{idx} = obj.Ribbons{idx}.CreateClone();
        end
        
        ret = CombineRibbons('Ribbons', Ribbons_cloned);
        
    end
    
    
    
%% getTotalHeight
%> @brief Gets the total height of the ribbon structure.
%> @return Returns with the total height (length) of the scattering region in units of the lattice vector.
    function ret = getTotalHeight( obj )
        heights = obj.getHeights();
        ret = sum( heights );
        
        if length(heights) == 1
            return
        end
        
        % accounting for the coupling between the ribbons
        if ~isempty( heights )
          ret = ret + length(obj.Ribbons )-1;
          end
        
        % accounting for the interface regions
        for idx = 2:length( obj.Ribbons )
            ribbon_tmp = obj.Ribbons{idx};
            if ~isempty(ribbon_tmp.Interface_Regions)
                ret = ret + 1;
            end
        end
    end    

       

%% addRibbon
%> @brief Adds a Ribbon part to the system.
%> @param ribbon2add An instance of the ribbon class.
%> @param varargin Cell array of optional parameters:.
function addRibbon( obj, ribbon2add, varargin )
    
    p = inputParser;
    p.addParameter('shiftcoordinates', true );   
    p.parse(varargin{:});       
    shiftcoordinates       = p.Results.shiftcoordinates;
    
    if ~isempty( obj.E )
          ribbon2add.setEnergy( obj.E );
    end
    
    % checking if the width of the ribbons is equal
    if obj.checkwidth( ribbon2add ) == -1
        err = MException('CombineRibbons:ribbon2add', 'widths of the ribbons are not identical.');
        save('Error_CombineRibbons_ribbon2add.mat');
        throw(err);
    end
    
    
    
    % shift the coordinates
    if shiftcoordinates
        total_height = obj.getTotalHeight();
        if ~isempty(obj.Ribbons)
          total_height = total_height + 1; % 1 for the coupling length
          if ~isempty(ribbon2add.Interface_Regions)
               total_height = total_height + 1;
          end
        end
        ribbon2add.ShiftCoordinates( total_height );
    end
    
    
    obj.Ribbons{end+1} = ribbon2add;    
   
    obj.height = obj.getTotalHeight();
    
end




end % methods public

methods ( Access = protected )




% getHeights
%> @brief Returns an array of the heights of the ribbon parts.
%> @return Returns Returns an array of the heights of the ribbon parts.
    function ret = getHeights( obj )
        ribbonnum = length( obj.Ribbons );
        ret = zeros( ribbonnum, 1); 
        for idx = 1:ribbonnum
            ribbon_tmp = obj.Ribbons{idx};
            ret(idx) =  ribbon_tmp.height;
        end
    end
    


    




    
%% attachRibbon
%> @brief Attach the next ribbon to the ribbon sequence for the calculation of the surface Greens function. 
%> @param ribbon_tmp An instance of the ribbon class.
%> @param varargin Cell array of optional parameters:.
function attachRibbon( obj, ribbon_tmp, varargin )
    
        p = inputParser;
        p.addParameter('onlyGinv', false );
        p.addParameter('decimate', true );
        p.addParameter('gfininvfromHamiltonian', false);
        p.addParameter('PotInScatter', []);
        
        p.parse(varargin{:});
       
        onlyGinv           = p.Results.onlyGinv;
        decimate           = p.Results.decimate;
        gfininvfromHamiltonian = p.Results.gfininvfromHamiltonian;
        PotInScatter = p.Results.PotInScatter;
        
            
    if isempty(obj.Ginv) && isempty(obj.G)
        if gfininvfromHamiltonian
            ribbon_tmp.CalcFiniteGreensFunctionFromHamiltonian( 'onlyGinv', true, 'gauge_trans', false, 'PotInScatter', PotInScatter );
        else
            ribbon_tmp.CalcFiniteGreensFunction( 'onlyGinv', true, 'gauge_trans', false );
        end
        [~,obj.Ginv] = ribbon_tmp.GetFiniteGreensFunction();
        obj.coordinates = ribbon_tmp.getCoordinates();
    else
        if isempty(obj.Ginv)
            if isnan(rcond(obj.G) ) || abs( rcond(obj.G) ) < 1e-15
                obj.Ginv = obj.inv_SVD( G );
            else
                obj.Ginv = inv(obj.G);
            end       
            obj.G = [];
        end
        size_gfin = size(obj.Ginv);
        if gfininvfromHamiltonian
            ribbon_tmp.CalcFiniteGreensFunctionFromHamiltonian( 'onlyGinv', true, 'gauge_trans', false );
        else
            ribbon_tmp.CalcFiniteGreensFunction( 'onlyGinv', true, 'gauge_trans', false );
        end
        [~,gfininv_tmp] = ribbon_tmp.GetFiniteGreensFunction();
        size_gfin_tmp = size( gfininv_tmp);
            
        Surface_sc = ribbon_tmp.Scatter_UC.CreateClone();
            
        if isempty( ribbon_tmp.Interface_Regions )
                
            H1_sc_tmp    = Surface_sc.Read( 'H1' );
            H1adj_sc_tmp = Surface_sc.Read( 'H1adj' );
            size_H1_sc_tmp = size(H1_sc_tmp);
                
            H1_sc = zeros(size_gfin(1), size_gfin_tmp(2) );
            H1adj_sc = zeros(size_gfin_tmp(1), size_gfin(2) );
                
            H1_sc(end-size_H1_sc_tmp(1)+1:end, 1:size_H1_sc_tmp(2)) = H1_sc_tmp;
            H1_sc_tmp = [];
                
            H1adj_sc( 1:size_H1_sc_tmp(2), end-size_H1_sc_tmp(1)+1:end) = H1adj_sc_tmp;
            H1adj_sc_tmp = [];
                
            obj.Ginv = [obj.Ginv, -H1_sc; -H1adj_sc, gfininv_tmp];
            
            coordinates_tmp = ribbon_tmp.getCoordinates();      
            fnames = fieldnames(obj.coordinates);
            for idx = 1:length( fnames )
                fname = fnames{idx};
                if strcmp( fname, 'a' )
                    continue
                else
                    obj.coordinates.(fname) = [obj.coordinates.(fname); coordinates_tmp.(fname)];
                end
            end
            
                
                
        else
                
            H0_interface = ribbon_tmp.Interface_Regions{1}.Read('H0');
            H1_interface_tmp = ribbon_tmp.Interface_Regions{1}.Read('H1');
            H1adj_interface_tmp = ribbon_tmp.Interface_Regions{1}.Read('H1adj');
            size_interface = size(H0_interface,1);
            size_H1_interface_tmp = size( H1_interface_tmp );
                
            % ribbon_tmp -> interface
            H1_sc_tmp    = Surface_sc.Read( 'H1' );
            H1adj_sc_tmp = Surface_sc.Read( 'H1adj' );
            size_H1_sc_tmp = size(H1_sc_tmp);
                
            H1_sc = zeros(size_interface, size_gfin_tmp(2) );
            H1adj_sc = zeros(size_gfin_tmp(1), size_interface );
                
            H1_sc(end-size_H1_sc_tmp(1)+1:end, 1:size_H1_sc_tmp(2)) = H1_sc_tmp;
            H1_sc_tmp = [];
                
            H1adj_sc( 1:size_H1_sc_tmp(2), end-size_H1_sc_tmp(1)+1:end ) = H1adj_sc_tmp;
            H1adj_sc_tmp = [];
                
            % interface -> other ribbons
            H1_interface = zeros(size_gfin(1), size_interface );
            H1adj_interface = zeros( size_interface, size_gfin(2) );
                
            H1_interface(end-size_H1_interface_tmp(1)+1:end, 1:size_H1_interface_tmp(2)) = H1_interface_tmp;
            H1_interface_tmp = []; 
                
            H1adj_interface(1:size_H1_interface_tmp(2), end-size_H1_interface_tmp(1)+1:end ) = H1adj_interface_tmp;
            H1adj_interface_tmp = [];                 
            
                
            obj.Ginv = [ obj.Ginv, -H1_interface, zeros(size_gfin(1),size_gfin_tmp(2)); ...
                  -H1adj_interface, obj.E*eye(size(H0_interface))-H0_interface, -H1_sc; ...
                zeros(size_gfin_tmp(1),size_gfin(2)), -H1adj_sc, gfininv_tmp]; 
                             

            [coordinates_tmp, coordinates_iterface] = ribbon_tmp.getCoordinates();
            fnames = fieldnames(obj.coordinates);
            for idx = 1:length( fnames )
                fname = fnames{idx};
                if strcmp( fname, 'a' )
                    continue
                else
                    obj.coordinates.(fname) = [obj.coordinates.(fname); coordinates_iterface{1}.(fname); coordinates_tmp.(fname)];
                end
            end
            
                
        end 
        
              
           
        if decimate
            size_gfininv = size(obj.Ginv,1);
            M = Surface_sc.Read('M');
            
            gfin_tmp = inv(obj.Ginv);
            gfin_tmp = [gfin_tmp(1:M, 1:M), gfin_tmp(1:M, size_gfininv-M+1:size_gfininv); ...
                                gfin_tmp(size_gfininv-M+1:size_gfininv ,1:M), gfin_tmp(size_gfininv-M+1:size_gfininv, size_gfininv-M+1:size_gfininv)];
            obj.Ginv = inv(gfin_tmp);
            
            names = fieldnames(obj.coordinates);
            for idx = 1:length(names)
                if strcmp(names{idx}, 'a') || strcmp(names{idx}, 'b')
                    continue;
                end
                
                array_tmp = obj.coordinates.(names{idx});
                if isempty(array_tmp)
                    continue;
                end
                array_tmp = array_tmp( [1:M, size_gfininv-M+1:size_gfininv]);
                obj.coordinates.(names{idx}) = array_tmp;
            end
            
%             kulso_szabfokok = [1:M, size_gfininv-M+1:size_gfininv];
%             [gfininv, coordinates] = ribbon.DecimationFunction( kulso_szabfokok, gfininv, 'coordinates', coordinates);
        end
            
        
    end
    
              
    
    if ~onlyGinv
        if isnan(rcond(obj.Ginv) ) || abs( rcond(obj.Ginv) ) < 1e-15
            obj.G = obj.inv_SVD( obj.Ginv );
        else
            obj.G = inv(obj.Ginv);
        end       
        obj.Ginv = [];
    else
        obj.G = [];
    end       
        
end




%% createInterfaceRegions
%> @brief Creates instances of class Surface_Green_function describing an interface region between the leads and the ribbon parts. The created interface is stored within the attribute Interface_Regions_all.
function createInterfaceRegions( obj )

    Idx = length(obj.Ribbons) + 1;
    obj.Interface_Regions_all = cell(Idx,1);
    
    % crating the interface region between the leads and the scattring region
    obj.Interface_Regions_all{1} = obj.Interface_Regions{1};
    obj.Interface_Regions_all{Idx} = obj.Interface_Regions{2};
    
    % creating the interface regions between the combined ribbons
    for idx = 1:length(obj.Ribbons)
        ribbon_tmp = obj.Ribbons{idx};
        if isempty( obj.Interface_Regions_all{idx+1} )
          obj.Interface_Regions_all{idx+1} = ribbon_tmp.Scatter_UC.CreateClone();
        end
        obj.Ribbons{idx}.setInterfaceRegions( {obj.Interface_Regions_all{idx}, obj.Interface_Regions_all{idx+1}} ); 
    end
    
    ribbon_tmp = obj.Ribbons{1};
    obj.Scatter_UC = ribbon_tmp.Scatter_UC;
    obj.CreateInterface(1);
    ribbon_tmp = obj.Ribbons{end};
    obj.Scatter_UC = ribbon_tmp.Scatter_UC;
    obj.CreateInterface(2);
    
    % Decimating the interface regions
    for idx = 2:length(obj.Ribbons)
        obj.Ribbons{idx}.CreateInterface(1)
    end  
   
        
end


%% checkwidth
%> @brief Checks the width of the added ribbon.
%> @param ribbon_tmp An instance of class Ribbon
%> @return Returns with true if the width of the added ribbon is consitent with the previous ones.
function ret = checkwidth( obj, ribbon_tmp )
    
    % checking if the width of the ribbons is equal
    if isempty(obj.width)
        obj.width = ribbon_tmp.width;
        ret = 0;
    elseif obj.width ~= ribbon_tmp.width
        ret = -1;
    else
        ret = 0;
    end
        
        
end


    end % protected methods
    
    
    methods (Access=private)    
   

%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Cell array of optional parameters:
%> @param 'Ribbons' Cell array of classes #Ribbon to be combined.
    function InputParsing( obj, varargin )
        
        p = inputParser;
          p.addParameter('Ribbons', []);
        
        p.parse(varargin{:});

          Ribbons      = p.Results.Ribbons;


        if ~isempty(Ribbons)
             for idx = 1:length(Ribbons)
                 ribbon_tmp = Ribbons{idx};
                 
                 if idx == 1
                     names = properties('Ribbon');
                     for jdx = 1:length(names);
                         obj.(names{jdx}) = ribbon_tmp.(names{jdx});
                     end
                     obj.Opt = ribbon_tmp.getOpt();
                     obj.param = ribbon_tmp.param;
                 else
                     % checking if the width of the ribbons is equal
                     if obj.checkwidth( ribbon_tmp ) == -1
                         err = MException('CombineRibbons:InputParsing', 'widths of the ribbons are not identical.');
                         save('Error_CombineRibbons_InputParsing.mat');
                         throw(err);
                     end
                 end
                 
                 
                 obj.addRibbon( ribbon_tmp );
            end
        else
            err = MException('CombineRibbons:InputParsing', 'obj.Ribbons must contain at least one Ribbon interface.');
            save('Error_CombineRibbons_InputParsing.mat');
            throw(err);            
        end
        
    end
    
end % private methods
end