IV.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - IV
%    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 IV.m
%> @brief A class to calculate the I-V curve for a two terminal setup, based on the non-equilibrium Greens function technique of Eur. Phys. J. B 53, 537-549 (2006)
%> @} 
%> @brief A class to calculate the I-V curve for a two terminal setup, based on the non-equilibrium Greens function technique of Eur. Phys. J. B 53, 537-549 (2006)
%> @Available
%> EQuUs v4.8 or later
%%
classdef IV < UtilsBase
  
    
    properties (Access = protected)
        %> The energy array used in the calculations
        Evec 
        %> The maximal bias to be used in the calculations
        bias_max 
        %> The number of the energy points in the attribute Evec
        Edb
    end
    
    
methods (Access=public)
        
%% Contructor of the class
%> @brief Constructor of the class.
%> @param Opt an instance of class #Opt.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
    function obj = IV( Opt, varargin )  
        
        obj = obj@UtilsBase(Opt, varargin{:} );

        if strcmpi(class(obj), 'IV') 
            obj.InputParsing( varargin{:});
        end
        
        
    end


    
%% CurrentCalc   
%> @brief Calculates the current-voltage relation as a function of the bias.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'tuned_contact' String value. Identifies the tuned contact for the calculations. Possible values are: 'right', 'left', 'symmetric'.
%> @return [1] The calculated current array.
%> @return [2] The calculated voltage bias array.
%> @return [3] The calculated differential conductance.
    function [current, bias, DifferentialConductance_vec] = IVcalc( obj, varargin )

          p = inputParser;
        p.addParameter('tuned_contact', 'right', @ischar); %which contact is tuned during the calculations: right, left, symmetric
        
        p.parse(varargin{:});
       
        tuned_contact = p.Results.tuned_contact;
        
        % create energy array if it is empty
        if isempty( obj.Evec )
            obj.CreateEnergyArray( 'tuned_contact', tuned_contact );
        end
        
            
        
        Evec_tmp = obj.Evec;
        
        % preallocating array for the differential conductance
        DifferentialConductance_vec = zeros( size(Evec_tmp) );
        
        % opening the parallel pool
        poolmanager = Parallel( obj.Opt );
        poolmanager.openPool(); 
        
        % 
        cNTerminal_loc = obj.junction;  
        
        % creating function handles for the parfor loop
        hDifferentialConductanceT0 = @(Energy, cNTerminal_loc)(obj.DifferentialConductanceT0(Energy, cNTerminal_loc));
        
        parfor (iidx = 1:length(Evec_tmp), poolmanager.NumWorkers)           
            Energy = Evec_tmp(iidx);                   
            DifferentialConductance_vec(iidx) = feval(hDifferentialConductanceT0, Energy, cNTerminal_loc);
        end
        
        % closing the parallel pool
        poolmanager.closePool();
        
        % Calculate the IV curve
        if strcmpi( tuned_contact, 'right') || strcmpi( tuned_contact, 'left')
            current = zeros( size(Evec_tmp) );
            for idx = 2:length(Evec_tmp)
                current(idx) = obj.currentcalc( DifferentialConductance_vec(1:idx), Evec_tmp(1:idx), tuned_contact );
            end    
            
            indexes = (0 <= Evec_tmp & Evec_tmp <= obj.bias_max);
            bias = Evec_tmp(indexes);
            current = current(indexes);
            DifferentialConductance_vec = DifferentialConductance_vec(indexes);            
                        
        elseif strcmpi( tuned_contact, 'symmetric' )
            Evec_length = floor(length( Evec_tmp )/2);
            current = zeros( Evec_length, 1 );
            for idx = 2:Evec_length
                current(idx) = obj.currentcalc( DifferentialConductance_vec(Evec_length-idx+1: Evec_length+idx-1), Evec_tmp(Evec_length-idx+1: Evec_length+idx-1), tuned_contact );
            end
            bias = [ 0, Evec_tmp(end-Evec_length+1:end) - Evec_tmp(Evec_length:-1:1)];    
            DifferentialConductance_vec = [DifferentialConductance_vec(Evec_length+1), DifferentialConductance_vec(Evec_length:-1:1) + DifferentialConductance_vec(end-Evec_length+1:end)];
            indexes = 1;
        end 
        
        
          


    end
    
        
    
    
end % public methods



methods (Access=protected)

%> @brief Calculates the current by integrating the differential conductance in a bias window
%> @param differentialConductance Array of differential conductance
%> @param Evec Array of energy values
%> @param tuned_contact String value. Identifies the tuned contact for the calculations. Possible values are: 'right', 'left', 'symmetric'.
%> @return Returns with the calculated current
     function current = currentcalc( obj, differentialConductance, Evec, tuned_contact )

          deltaEvec = [0, diff(Evec)];
        
        % defining the Fermi functions for the individual leads
        if strcmpi( tuned_contact, 'right') 
            Fermi_left  = @(E)(obj.Fermi( E ));
            Fermi_right = @(E)(obj.Fermi( E-obj.bias_max ));
        elseif strcmpi( tuned_contact, 'left') 
            Fermi_left  = @(E)(obj.Fermi( E-obj.bias_max ));
            Fermi_right = @(E)(obj.Fermi( E ));
        elseif strcmpi( tuned_contact, 'symmetric') 
            Fermi_left  = @(E)(obj.Fermi( E+obj.bias_max/2 ));
            Fermi_right = @(E)(obj.Fermi( E-obj.bias_max/2 ));            
        else
            current = [];
            return
        end
        
        % including temperature
        differentialConductance = differentialConductance.*(Fermi_right(Evec) - Fermi_left(Evec));
        
        if mod(length( obj.Evec ),2) == 0
            current = obj.SimphsonInt( differentialConductance(1:end-1).*deltaEvec(1:end-1) );
            current = current + (differentialConductance(end-1) + differentialConductance(end) )*deltaEvec(end)/2;
        else
            current = obj.SimphsonInt( differentialConductance.*deltaEvec );            
        end 

     end
    
    
    
%% create_scatter_GreensFunction
%> @brief Calculates the surface Green operator of the scattering region.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'gauge_trans' Logical value. Set true (default) to perform gauge transformation on the Green's function and on the Hamiltonians, or false otherwise.
%> @param 'junction' An instance of class #NTerminal or its subclass.
%> @return The inverse of the Green operator.
    function gfininv = create_scatter_GreensFunction(  obj, varargin )
        
        p = inputParser;
        p.addParameter('gauge_trans', true); % logical: true if want to perform gauge transformation on the Green's function and Hamiltonians
        p.addParameter('junction', obj.junction); %for parallel computations
        p.parse(varargin{:});       
        gauge_trans     = p.Results.gauge_trans;
        junction_loc     = p.Results.junction;       
        
        if obj.gfininvfromHamiltonian
                %calculating the Greens function from the Hamiltonian
                junction_loc.CalcFiniteGreensFunctionFromHamiltonian('onlyGinv', true, 'gauge_trans', false, 'scatterPotential', obj.scatterPotential );
        else
                % Calculating the Greens function by the fast way optimized for translational invariant systems
                junction_loc.CalcFiniteGreensFunction('onlyGinv', true, 'gauge_trans', false);
        end
        [~, gfininv] = junction_loc.GetFiniteGreensFunction();
        coordinates_scatter = junction_loc.getCoordinates();
               
        
        % gauge transformation of the vector potential in the effective Hamiltonians
        if gauge_trans && ~isempty(  junction_loc.PeierlsTransform_Scatter )
            gfininv = junction_loc.PeierlsTransform_Scatter.gaugeTransformation( gfininv, coordinates_scatter, junction_loc.gauge_field );
            %NTerminal_loc.PeierlsTransform_Scatter.gaugeTransformationOnLead( NTerminal_loc.Surface_tmp, NTerminal_loc.gauge_field );  
        end
            
        
    end
       
    
end
    
    
methods (Access=protected)  %NEW

%% DifferentialConductanceT0
%> @brief Calculates the zero temperature differential conductance for a given energy
%> @param Energy The energy value (use the same units as in the Hamiltonian)
%> @param junction An instance of class #NTerminal or its subclass.
%> @return Returns with the calculated zero temperature differential conductance.
     function ret = DifferentialConductanceT0( obj, Energy, junction )
        
            % setting the energy value
          obj.SetEnergy( Energy, 'junction', junction );
               cTransport_Interface = junction.FL_handles;
        
            % calculating the Green function of the scattering reigon
          gfininv = obj.create_scatter_GreensFunction('junction', junction);         
        
            % setting the function handle for the Dyson equation
          Dysonfunc = @()(junction.CustomDysonFunc( 'decimate', true, 'gfininv', gfininv, 'constant_channels', false, ...
                'SelfEnergy', obj.useSelfEnergy) );

               cTransport_Interface.DysonEq( 'CustomDyson', Dysonfunc );
       
          
            ret = cTransport_Interface.Conductance2();

    end
    
    
    
%% CreateEnergyArray
%> @brief Creates the array of energy values for the integration
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'tuned_contact' String value. Identifies the tuned contact for the calculations. Possible values are: 'right', 'left', 'symmetric'.
     function CreateEnergyArray( obj, varargin ) 

        p = inputParser;
        p.addParameter('tuned_contact', 'right', @ischar); %which contact is tuned during the calculations: right, left, symmetric
        
        p.parse(varargin{:});
       
        tuned_contact = p.Results.tuned_contact;
        
        if ~isscalar( obj.T )
            error( 'EQuUs:Utils:IV:CreateEnergyArray', 'Temperature should be a scalar.');
        end
        
        if obj.T < obj.T_treshold
        
            if strcmpi( tuned_contact, 'right') 
                mu_left = 0;
                mu_right = obj.bias_max;
                Emin = min( [mu_left, mu_right]);
                Emax = max( [mu_left, mu_right]);
                obj.Evec = Emin:(Emax-Emin)/(obj.Edb+1):Emax;
                return
                
            elseif strcmpi( tuned_contact, 'left')
                mu_left = obj.bias_max;
                mu_right = 0;                
                Emin = min( [mu_left, mu_right]);
                Emax = max( [mu_left, mu_right]);
                obj.Evec = Emin:(Emax-Emin)/(obj.Edb+1):Emax;
                return
                
            elseif strcmpi( tuned_contact, 'symmetric' )
                mu_left = -obj.bias_max/2;
                mu_right = obj.bias_max/2;
                Emin = min( [mu_left, mu_right]);
                Emax = max( [mu_left, mu_right]);
                Emean = (Emax+Emin)/2;
                deltaEvec = (Emax-Emin)/obj.Edb;
                obj.Evec = [sort(Emean-deltaEvec:-deltaEvec:Emin, 'ascend'), Emean:deltaEvec:Emax];
                return
                
            else
                error( 'EQuUs:Utils:IV:CreateEnergyArray', 'Undefined contact to be tuned');
            end
            
        else
            tolerance = 1e-5;            
            
            if strcmpi( tuned_contact, 'right')
                mu_left = 0;
                mu_right = obj.bias_max;    
                Fermi_left  = @(E)(obj.Fermi( E-mu_left ));
                Fermi_right = @(E)(obj.Fermi( E-mu_right ));
                Emin = fzero(@(x)(Fermi_left(x)-1+tolerance), mu_left);
                Emax = fzero(@(x)(Fermi_right(x)-tolerance), mu_right);                
                obj.Evec = Emin:(Emax-Emin)/(obj.Edb+1):Emax;
                return
                
            elseif strcmpi( tuned_contact, 'left')
            
                mu_left = obj.bias_max; %relative to the Fermi level
                mu_right = 0;   
                Fermi_left  = @(E)(obj.Fermi( E-mu_left ));
                Fermi_right = @(E)(obj.Fermi( E-mu_right ));                
                Emin = fzero(@(x)(Fermi_left(x)-tolerance), mu_left);
                Emax = fzero(@(x)(Fermi_right(x)-1+tolerance), mu_right);    
                obj.Evec = Emin:(Emax-Emin)/(obj.Edb+1):Emax;
                return 
                
            elseif strcmpi( tuned_contact, 'symmetric' )
                mu_left = obj.junction.getFermiEnergy() -obj.bias_max/2;
                mu_right = obj.junction.getFermiEnergy() + obj.bias_max/2;
                Fermi_left  = @(E)(obj.Fermi( E-mu_left));
                Fermi_right = @(E)(obj.Fermi( E-mu_right ));                
                if obj.bias_max > 0
                    Emin = fzero(@(x)(Fermi_left(x)-1+tolerance), mu_left);
                    Emax = fzero(@(x)(Fermi_right(x)-tolerance), mu_right);
                elseif obj.bias_max < 0
                    Emin = fzero(@(x)(Fermi_left(x)-tolerance), mu_left);
                    Emax = fzero(@(x)(Fermi_right(x)-1+tolerance), mu_right);    
                else
                    obj.Evec = [];
                    return
                end
                
                Emean = (Emax+Emin)/2;
                deltaEvec = (Emax-Emin)/obj.Edb;
                obj.Evec = [sort(Emean-deltaEvec:-deltaEvec:Emin, 'ascend'), Emean:deltaEvec:Emax];
                return
            else
                error( 'EQuUs:Utils:IV:CreateEnergyArray', 'Undefined contact to be tuned');
            end 
            
        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 'T' The temperature in Kelvin (scalar or an array) 
%> @param 'scatterPotential' A function handle pot=f( #coordinates ) or pot=f( #CreateHamiltonians, Energy) for the potential to be applied in the Hamiltonian (used when FiniteGreensFunctionFromHamiltonian=true).
%> @param 'gfininvfromHamiltonian' Logical value. Set true calculate the surface Greens function of the scattering region from the Hamiltonaian of the scattering region, or false (default) to calculate it by the fast way (see Phys. Rev. B 90, 125428 (2014) for details).
%> @param 'useSelfEnergy' 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 'junction' An instance of class #NTerminal (or its derived class) representing the junction.
%> @param 'bias_max' The maximal bias to be used in the calculations
%> @param 'Edb' The number of the energy points in the attribute #Evec 
    function InputParsing( obj, varargin )
        
        
        p = inputParser;
        p.addParameter('T', obj.T, @isscalar);
        p.addParameter('junction', obj.junction); %NEW changed parameter name
        p.addParameter('scatterPotential', []);
        p.addParameter('gfininvfromHamiltonian', false);   
        p.addParameter('useSelfEnergy', true);    
        p.addParameter('bias_max', 0);  % The maximal bias to be used in the calculations %NEW
        p.addParameter('Edb', 300);  % The number of the energy points in the attribute Evec        %NEW
        
        p.parse(varargin{:});
        
        InputParsing@UtilsBase( obj, 'T', p.Results.T, ...
                            'junction', p.Results.junction, ...
                            'scatterPotential', p.Results.scatterPotential, ...
                            'gfininvfromHamiltonian', p.Results.gfininvfromHamiltonian, ...
                            'useSelfEnergy', p.Results.useSelfEnergy);
         
        obj.bias_max            = p.Results.bias_max;
        obj.Edb                 = p.Results.Edb;       
        
        if obj.bias_max <= 0
            error(['EQuUs:Utils:', class(obj.junction), ':InputParsing'], 'Maximal bias should be greater than zero.');
        end
               

    end        
        
        
end % protected methods    
    
end