NTerminal_Keldysh.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - NTerminal_Keldysh
%    Copyright (C) 2016 Peter Rakyta, Ph.D.
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see http://www.gnu.org/licenses/.
%
%> @addtogroup utilities Utilities
%> @{
%> @file NTerminal_Keldysh.m
%> @brief  A class describing an N-terminal geometry for steady state non-equilibrium calculations.
%> @image html two_terminal_structure.jpg
%> @image latex two_terminal_structure.jpg
%> @} 
%> @brief  A class describing an N-terminal geometry for steady state non-equilibrium calculations.
%> @Available
%> EQuUs v4.9 or later
%> <tr class="heading"><td colspan="2"><h2 class="groupheader"><a name="avail"></a>Structure described by the class</h2></td></tr>
%> @image html two_terminal_structure.jpg
%> @image latex two_terminal_structure.jpg
%> The drawing represents a two-terminal structure made of two leads, of a scattering region and of two interface regions between the leads and scattering center.
%> However, additional lead can be added to the structure.
%> Each rectangle describes a unit cell including singular and non-singular sites. 
%> The scattering center is, on the other hand, represented by a larger rectangle corresponding to the Hamiltonian of the scattering region.
%> Arrows indicate the hopping direction stored by the attributes in the corresponding classes (see attributes #CreateLeadHamiltonians.H1 and #InterfaceRegion.Hcoupling for details).
%> The orientation of the lead is +1 if the lead is terminated by the interface region in the positive direction, and -1 if the lead is terminated by the interface region in the negative direction.
%> (see attribute #CreateLeadHamiltonians.Lead_Orientation for details)
%%
classdef NTerminal_Keldysh < NTerminal & FermiDirac

    properties (Access = protected)
        %> An array containing the bias of the leads in the same unit as other energy scales in the Hamiltonians.
        bias_leads
        %> Lesser Green operator projected onto the scattering region
        lesserG 
    end       
    
    properties (Access = public)
        
    end
    


methods ( Access = public )
    
    
%% Contructor of the class
%> @brief Constructor of the class.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
function obj = NTerminal_Keldysh( varargin )
    obj = obj@NTerminal();
    obj = obj@FermiDirac();
    
    
    obj.lesserG = [];
    obj.bias_leads = [];
    
    
    if strcmpi( class(obj), 'NTerminal_Keldysh') 
        obj.InputParsing( varargin{:} );
    
        obj.cCustom_Hamiltonians = [];

        obj.display(['EQuUs:Utils:', class(obj), ':Creating NTerminal_Keldysh class.'])
    
        % exporting calculation parameters into an XML format
        createOutput( obj.filenameOut, obj.Opt, obj.param );            
            
        % initializing attributes
        obj.CreateHandles();
    
        % create Hamiltonian of the scattering region (or load from external source)
        obj.CreateNTerminalHamiltonians();
    end
    


end


%% CustomDysonFunc
%> @brief Custom Dyson function for a two terminal arrangement on a two dimensional lattice.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'gfininv' The inverse of the Greens function of the scattering region. For default the inverse of the attribute #G is used.
%> @param 'constant_channels' Logical value. Set true (default) to keep constant the number of the open channels in the leads for each energy value, or false otherwise.
%> @param 'onlyGinverz' Logical value. Set true to calculate only the inverse of the total Green operator, or false (default) to calculate #G as well.
%> @param 'recalculateSurface' A vector of the identification numbers of the lead surfaces to be recalculated.
%> @param 'decimate' Logical value. Set true (default) to eliminate all inner sites in the Greens function and keep only the selected sites. Set false to omit the decimation procedure.
%> @param 'kulso_szabfokok' Array of sites to be kept after the decimation procedure. (Use parameter 'keep_sites' instead)
%> @param 'selfEnergy' Logical value. Set true to use the self energies of the leads in the Dyson equation, or false (default) to use the surface Green function instead.
%> @param 'keep_sites' Name of sites to be kept in the resulted Green function (Possible values are: 'scatter', 'interface', 'lead').
%> @return [1] The calculated Greens function. 
%> @return [2] The inverse of the Green operator. 
%> @return [3] An instance of structure #junction_sites describing the sites in the calculated Green operator.
    function [Gret, Ginverz, junction_sites] = CustomDysonFunc( obj, varargin ) %NEW output
        
    p = inputParser;
    p.addParameter('gfininv', []);
    p.addParameter('constant_channels', true);
    p.addParameter('onlyGinverz', false );
    p.addParameter('recalculateSurface', 1:length(obj.param.Leads) );
    p.addParameter('decimate', true );
    p.addParameter('kulso_szabfokok', []); %The list of sites to be left after the decimation procedure
    p.addParameter('SelfEnergy', false); %set true to calculate the Dyson equation with the self energy
    p.addParameter('keep_sites', 'lead'); %Name of sites to be kept (scatter, interface, lead)
    
    p.parse(varargin{:});
    gfininv     = p.Results.gfininv;
    constant_channels = p.Results.constant_channels;
    onlyGinverz        = p.Results.onlyGinverz;
    recalculateSurface = p.Results.recalculateSurface;  
    decimate           = p.Results.decimate; 
    kulso_szabfokok = p.Results.kulso_szabfokok;
    useSelfEnergy         = p.Results.SelfEnergy;
    keep_sites        = p.Results.keep_sites;
    

    if ~isempty(recalculateSurface)
    
        % creating interfaces for the Leads
        if constant_channels
            shiftLeads = ones(length(obj.param.Leads),1)*obj.E;
        else
            shiftLeads = ones(length(obj.param.Leads),1)*0;
        end
    
        % creating Lead instaces and calculating the retarded surface Green operator/self-energy
        obj.FL_handles.LeadCalc('shiftLeads', shiftLeads, ...
            'SelfEnergy', useSelfEnergy, 'SurfaceGreensFunction', ~useSelfEnergy, 'gauge_field', obj.gauge_field, 'leads', recalculateSurface, 'q', obj.q, ...
            'leadmodel', obj.leadmodel, 'bias_leads', obj.bias_leads, 'T', obj.T, 'CustomHamiltonian', obj.cCustom_Hamiltonians);
        
        for idx = 1:length(recalculateSurface)
            obj.CreateInterface( recalculateSurface(idx) );
        end
        
    end
    
    [Gret, Ginverz, junction_sites] = CustomDysonFunc@NTerminal(obj, 'gfininv', gfininv, ...
                                    'onlyGinverz', onlyGinverz, ...
                                    'recalculateSurface', [], ...
                                    'decimate', decimate, ...
                                    'kulso_szabfokok', kulso_szabfokok, ...
                                    'SelfEnergy', useSelfEnergy, ...
                                    'keep_sites', keep_sites);                                 

        
    end

%% LesserGreenFunction
%> @brief Calculates the lesser Green function projected to the surface sites of the central device and the interface regions for eregy given #E. 
%> (Surface sites of the central device are directly connected to the interface regions.) The calculated Green function is stored by attributes #lesserG and #lesserGinv.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'gfininv' The inverse of the Greens function of the scattering region. For default the inverse of the attribute #G is used.
%> @param 'constant_channels' Logical value. Set true (default) to keep constant the number of the open channels in the leads for each energy value, or false otherwise.
%> @param 'recalculateSurface' A vector of the identification numbers of the lead surfaces to be recalculated.
%> @return [1] The calculated lesser Greens operator. 
%> @return [2] An instance of structure #junction_sites describing the sites in the calculated Green operator.
    function [Gret, junction_sites] = LesserGreenFunction( obj, varargin ) 
        
        p = inputParser;
        p.addParameter('gfininv', []);
        p.addParameter('constant_channels', false);
        p.addParameter('recalculateSurface', 1:length(obj.param.Leads) );
        p.parse(varargin{:});
        gfininv     = p.Results.gfininv;
        constant_channels = p.Results.constant_channels;
        recalculateSurface = p.Results.recalculateSurface; 
        
        % evaluate the Dyson equation for the retarded Green operator
        Dysonfunc = @()(obj.CustomDysonFunc( ...
                'decimate', false, 'constant_channels', constant_channels, ...
                'gfininv', gfininv, 'recalculateSurface', recalculateSurface, ...
                'SelfEnergy', true));
        [retardedG, junction_sites] = obj.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
        
        % create a global matrix containing the lesser self-energies
        % and obtaining the lesser self energies of the leads
        lesser_self_energy = zeros( size(retardedG) );
        Leads = obj.FL_handles.Read( 'Leads' );
        leadnum = length( Leads );
        for idx = 1:leadnum
            indexes = junction_sites.Leads{idx}.site_indexes;
            lesser_self_energy( indexes, indexes ) = Leads{idx}.LesserSelfEnergy();
        end
        
        % Calculate the lesser Green operator via
        % Eq (61) in Adv. Nat. Sci.: Nanosci. Nanotechnol. 5 (2014) 033001
        Gret = retardedG*lesser_self_energy*retardedG';
        
        
    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@NTerminal( obj, Energy );

        % resetting the class describing the N-terminal geometry
        if ~isempty( obj.FL_handles ) && strcmpi(class(obj.FL_handles), 'Transport_Interface_Keldysh')        
            obj.FL_handles.setEnergy( obj.E ); 
        end                 
        
        % recreate the Hamiltonian of the scattering region
        if ~isempty( obj.CreateH ) && strcmpi( class(obj), 'NTerminal_Keldysh' )      
            obj.CreateNTerminalHamiltonians();
        end
        
    end    
    
    
%% setTemperature
%> @brief Sets the temperature in the calculations
%> @param T The value of the temperature in Kelvins.
    function setTemperature( obj, T )
        setTemperature@FermiDirac( obj, T );

        % setting the temperature in class describing the N-terminal geometry
        if ~isempty( obj.FL_handles ) && strcmpi(class(obj.FL_handles), 'Transport_Interface_Keldysh')        
            obj.FL_handles.setTemperature( obj.T ); 
        end  
        
        obj.lesserG = [];
        
    end        


end % end methods public



methods (Access=protected)   
    
%% CreateHandles
%> @brief Initializes the attributes of the class.
    function CreateHandles( obj )  
        
        % creating classes for managing the magnetic field
        obj.setHandlesForMagneticField();
        
        % create class storing the Hamiltonian of the scattering region
        obj.CreateH = CreateHamiltonians_Keldysh(obj.Opt, obj.param, 'q', obj.q);
        
        % create class for transport calculations specific for steady state non-equilibrium calculations
        obj.FL_handles = Transport_Interface_Keldysh(obj.E, obj.Opt, obj.param, 'CreateH', obj.CreateH, 'PeierlsTransform', obj.PeierlsTransform_Leads );
        
        % read out structure Opt containing the computational parameters
        obj.Opt = obj.FL_handles.Read( 'Opt' );
        
        % creating classes of the interface regions
        obj.Interface_Regions = cell(length(obj.param.Leads),1);
        for idx = 1:length(obj.Interface_Regions)
            Lead_Orientation = obj.param.Leads{idx}.Lead_Orientation;
            obj.Interface_Regions{idx} = InterfaceRegion(obj.Opt, obj.param, 'Hanyadik_Lead', idx, 'q', obj.q, 'Lead_Orientation', Lead_Orientation);
        end
        
    end


%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'filenameIn' The input filename containing the computational parameters. (Use parameters 'Op' and 'param' instead)
%> @param 'filenameOut' The output filename to export the computational parameters.
%> @param 'WorkingDir' The absolute path to the working directoy.
%> @param 'CustomHamiltoniansHandle' function handle for the custom Hamiltonians. Has the same inputs as #Custom_Hamiltonians.LoadHamiltonians and output values defined by the example #Hamiltonians.
%> @param 'E' The energy value used in the calculations (in the same units as the Hamiltonian).
%> @param 'EF' The Fermi energy in the same units as the Hamiltonian. Attribute #E is measured from this value. (Use for equilibrium calculations in the zero temperature limit. Overrides the one comming from the external source)
%> @param 'silent' Set true to suppress output messages.
%> @param 'leadmodel' A function handle #Lead=f( idx, E, varargin ) of the alternative lead model with equivalent inputs and return values as #Transport_Interface.SurfaceGreenFunctionCalculator and with E standing for the energy.
%> @param 'interfacemodel' A function handle f( #InterfaceRegion ) to manually adjus the interface regions. (Usefull when 'leadmodel' is also given. For example see @InterfaceModel)
%> @param 'Opt' An instance of the structure #Opt.
%> @param 'param' An instance of the structure #param.
%> @param 'q' The transverse momentum quantum number.
%> @param 'T' The absolute temperature in Kelvins.
%> @param 'mu' The Chemical potential in the same unit as other energy scales in the Hamiltonians.
%> @param 'mu_leads' An array containing the chemical potentials of the leads in the same unit as other energy scales in the Hamiltonians.
    function InputParsing( obj, varargin)
    
        p = inputParser;
        p.addParameter('filenameIn', obj.filenameIn, @ischar);
        p.addParameter('filenameOut', obj.filenameOut, @ischar);
        p.addParameter('WorkingDir', obj.WorkingDir, @ischar);
        p.addParameter('CustomHamiltoniansHandle', obj.CustomHamiltoniansHandle); %function handle for the custom Hamiltonians
        p.addParameter('E', obj.E, @isscalar);
        p.addParameter('silent', obj.silent);   
        p.addParameter('leadmodel', obj.leadmodel); %individual physical model for the contacts
        p.addParameter('interfacemodel', obj.interfacemodel); %individual physical model for the interface regions
        p.addParameter('Opt', obj.Opt);
        p.addParameter('param', obj.param);
        p.addParameter('q', obj.q);
        p.addParameter('T', obj.T);
        p.addParameter('EF', obj.EF);
        p.addParameter('bias_leads', obj.bias_leads);
        
        p.parse(varargin{:});
        
        
        InputParsing@NTerminal( obj, 'filenameIn', p.Results.filenameIn, ...
            'filenameOut', p.Results.filenameOut, ...
            'WorkingDir', p.Results.WorkingDir, ...
            'CustomHamiltoniansHandle', p.Results.CustomHamiltoniansHandle, ...
            'E', p.Results.E, ...
            'EF', p.Results.EF, ...
            'silent', p.Results.silent, ...
            'leadmodel', p.Results.leadmodel, ...
            'interfacemodel', p.Results.interfacemodel, ...
            'q', p.Results.q, ...
            'Opt', p.Results.Opt, ...
            'param', p.Results.param);
        
        InputParsing@FermiDirac( obj, 'T', p.Results.T, ...
                    'mu', 0);
                
        
        % setting the bias in the leads
        if isempty( p.Results.bias_leads )
            obj.bias_leads = zeros( size( obj.param.Leads ) );
        else
            obj.bias_leads = p.Results.bias_leads;
        end
        
    end

end % methods private

end