QuantumDot.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - QuantumDot
%    Copyright (C) 2009-2017 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 QuantumDot.m
%> @brief A class to process transport calculations on quantum dots. Based on the real-time diagrammatic method of PRB 50, 18439
%> @} 
%> @brief A class to process transport calculations on quantum dots. Based on the real-time diagrammatic method of PRB 50, 18439
%> @Available
%> EQuUs v4.9 or later
%%
classdef QuantumDot < Density & IV
    
    properties ( Access = public )
        %> The charging energy
        Ec
    end    
    
    
    properties ( Access = protected )      
        %> Array of partition functions of 1<N<StatesNum occupied one-electron states ( Z = sum_n(exp(-beta*En)) )
        PartitionFunctionsOneElectron
        %> Occupation probabilities of an isolated, the N-occupied QD in thermal equilibrium
        ThermalOccupationProbabilities
        %> Occupation probabilities of the QD connected to the leads occupied by N-N0 excess charge
        OccupationProbabilities
        %> The calculated transition rates between the individual N-occupations states of the QD
        SelfEnergy
        %> The calculated tunneling rates between the individual leads and the QD
        TunnelingRates
        %> The calculated current operator between the QD and a chosen lead. (see method DetermineCurrentOperator)
        CurrentOp
        %> Array containing the occupation  numbers
        OccNumbers
        %> An array containing the chamical potentials corresponding to the occupation numbers of OccNumbers
        mu_QD
        %> Number of states utilized in the calculations.
        StatesNum
        
    end
    
    
    methods (Access=public)
    
%% QuantumDot
%> @brief Constructor of the class.
%> @param Opt An instance of structure #Opt.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
    function obj = QuantumDot( Opt, varargin )    
        
        
        obj = obj@Density( Opt, varargin{:} );  
        obj = obj@IV( Opt, varargin{:} );
        
        
        obj.Initialize();
        
        
        
    end
    
%% DetermineDensityOfStates
%> @brief Calculates the DOS to deremine the one-electron occupation levels. This must be done manually.
%> @param Emin The minimum of the energy array.
%> @param Emax The maximum of the energy array.
%> @param Edb The number of the energy points in the energy array.
    function DetermineDensityOfStates( obj, Emin, Emax, Edb )  
        
        %creating the energy array
        Evec = Emin:(Emax-Emin)/Edb:Emax;
        
        % calculating the desnity of states
        obj.DOSCalc(Evec);
        
        
           figure
           plot( obj.DOSdata.Evec, obj.DOSdata.DOS )
    end  
    

    
%% DetermineThermalOccupationProbabilities
%> @brief Determines the equilibrium occupation probabilities PN (see Eq (8) in Eur. Phys. J B 10, 119)
%> @param gateVoltage The back gate voltage in units of the energy unit used in the Hamiltonian parameters.
    function DetermineThermalOccupationProbabilities( obj, gateVoltage )  
        
        % determine partition functions for the occupied one-electron states
        % obj.DeterminePartitionFunction();
        
        N0 = length(obj.OccNumbers);
        
        % sum for the occupation numbers        
        OccupationProbabilities_tmp = zeros( N0, 1);
        for idx = 1:N0       
            total_charging_energy = obj.Ec/2*(obj.OccNumbers(idx) - gateVoltage/obj.Ec)^2;
            
            %OccupationProbabilities(idx) = obj.PartitionFunctionsOneElectron(idx) * exp( -obj.beta*total_charging_energy );
            OccupationProbabilities_tmp(idx) = exp( -obj.beta*total_charging_energy );
        
        end
        
        OccupationProbabilities_tmp = OccupationProbabilities_tmp/sum(OccupationProbabilities_tmp);        
        obj.ThermalOccupationProbabilities = OccupationProbabilities_tmp;

        
        
    end  
    
%% DetermineOccupationProbabilities
%> @brief Determines the non-equilibrium occupation probabilities P_N (see Eq (8) in Eur. Phys. J B 10, 119)
%> @return Returns with the calculated occupation probabilities
    function ret = DetermineOccupationProbabilities( obj )   
        
        SelfEnergy_tmp = obj.SelfEnergy;
        
        % an arbitrary chosen index. See Eq (6) in PRB 68, 115105
        idx0 = 1;
        
        SelfEnergy_tmp(:,idx0) = 1;
        
        obj.OccupationProbabilities = inv(SelfEnergy_tmp);
        obj.OccupationProbabilities = obj.OccupationProbabilities(idx0, :);
        
        % check the definition of the probability distribution
        tolerance = 1e-6;
        if sum( obj.OccupationProbabilities < -tolerance ) > 0
            error('EQuUs:Utils:QuantumDot:DetermineOccupationProbabilities', 'Negative probability obtained during the calculations.')
        end
        
        if abs(sum(obj.OccupationProbabilities) - 1) > tolerance
            error('EQuUs:Utils:QuantumDot:DetermineOccupationProbabilities', 'The sum of the probabilities is not equal to one.')
        end

          ret = obj.OccupationProbabilities;
        
%        disp(obj.OccupationProbabilities)

        
    end    
    
%% DetermineCurrent
%> @brief Calculates the current (in units of e^2/hbar) through a lead leadnum at gate voltage gateVoltage and at bias bias.
%> @param bias The bias voltage applied through the leads in units of the energy unit used in the Hamiltonian parameters.
%> @param gateVoltage The back gate voltage in units of the energy unit used in the Hamiltonian parameters.
%> @param leadnum Identification number of the lead. (see attribute #CreateLeadHamiltonians.Hanyadik_Lead).
%> @return Returns with the calculated current in units of \f$ e^2/\hbar \f$
    function Current = DetermineCurrent( obj, bias, gateVoltage, leadnum ) 
        
        % calculating the tunneling current
        if isempty( obj.CurrentOp )
            obj.DetermineCurrentOperator( bias, gateVoltage, leadnum );
        end 
        
        % calculating the occupation probabilities
        if isempty( obj.OccupationProbabilities )
            obj.DetermineOccupationProbabilities();
        end         
        
        % calculating the current
        Current = sum( obj.OccupationProbabilities*obj.CurrentOp ); 
    
        
        % Check the imaginary part of the current
        tolerance = 1e-8;
        if abs(Current) > tolerance && abs( imag(Current)/real(Current)) > tolerance
            warning('EQuUs:Utils:QuantumDot:DetermineCurrent', 'The current has to large imaginary part')
        end
        
        Current = real(Current);
        
    end
    
%% DetermineCurrentOperator
%> @brief Determine the current operator (in units of \f$ e^2/\hbar \f$) to calculate the current through a lead leadnum at gate voltage gateVoltage and at bias bias.
%> @param bias The bias voltage applied through the leads in units of the energy unit used in the Hamiltonian parameters.
%> @param gateVoltage The back gate voltage in units of the energy unit used in the Hamiltonian parameters.
%> @param leadnum Identification number of the lead. (see attribute #CreateLeadHamiltonians.Hanyadik_Lead).
    function DetermineCurrentOperator( obj, bias, gateVoltage, leadnum ) 

        bias_lead = (-1)^leadnum*bias/2;
        
        % calculating the tunneling rates
        if isempty( obj.TunnelingRates )
            obj.CalcTunnelingRates( gateVoltage );
        end

        
        % calculating Fermi-Dirac distribution functions in the chosen lead
        FD_lead_p = obj.Fermi(obj.Evec - bias_lead);
        FD_lead_m = 1 - FD_lead_p;

        % preallocating array
        CurrentOp = zeros( obj.StatesNum, obj.StatesNum);
        
        
        % calculate the transition rates (only the first order of the perturbation series)
        % from jdx to idx
        % The first order matrix elements are calculated by the diagrammatic rules in PRB 50, 18439
        deltaEvec = [0, diff(obj.Evec)]; 
        for idx = 1:obj.StatesNum             
            
            for jdx = 1:obj.StatesNum
                if idx-1==jdx        
                    
                    % calculating Fermi-Dirac distribution functions in the QD
                    mu_N = obj.mu_QD(jdx);
                    FD_QD_p = obj.Fermi(obj.Evec-mu_N-obj.Ec);
                    FD_QD_m = 1-FD_QD_p;      
                    
                    %Calculating the matrix elements
                    tmp = - 2*1i*obj.TunnelingRates{leadnum}(idx,:).*FD_lead_p.*FD_QD_m;
                    tmp = obj.SimphsonInt( tmp.*deltaEvec );
                    CurrentOp(jdx, idx) = CurrentOp(jdx, idx) + tmp;
                    
                elseif idx+1 == jdx
                    
                    % calculating Fermi-Dirac distribution functions in the QD
                    mu_N = obj.mu_QD(jdx);
                    FD_QD_p = obj.Fermi(obj.Evec-mu_N);
                    FD_QD_m = 1-FD_QD_p;       
                    
                    %Calculating the matrix elements
                    tmp = 2*1i*obj.TunnelingRates{leadnum}(jdx,:).*FD_lead_m.*FD_QD_p;
                    tmp = obj.SimphsonInt( tmp.*deltaEvec );
                    CurrentOp(jdx, idx) = CurrentOp(jdx, idx) + tmp;
                    
                else
                   continue
                end               
            end
        end
        
        CurrentOp = -1i*CurrentOp; % in units of e^2/hbar
        
        obj.CurrentOp = CurrentOp;
        
          
        
    end
    
%% DetermineSelfEnergy
%> @brief Determine the self-energy matrix elements describing the transition rates between the N occupated states of the QD. (see the diagrammatic rules in PRB 50, 18439)
%> @param bias The bias voltage applied through the leads in units of the energy unit used in the Hamiltonian parameters.
%> @param gateVoltage The back gate voltage in units of the energy unit used in the Hamiltonian parameters.
%> @return The calulated self-energy (transition rates between the N occupated states of the QD.)
    function SelfEnergy = DetermineSelfEnergy( obj, bias, gateVoltage ) 
        
        bias_lead = [ -bias/2  bias/2];
        
        % calculating the tunneling rates
        if isempty( obj.TunnelingRates )
            obj.CalcTunnelingRates( gateVoltage );
        end
        

        
        % calculating Fermi-Dirac distribution functions in the leads
        FD_lead_p = cell( length( obj.junction.Interface_Regions ), 1);
        FD_lead_m = cell( length( obj.junction.Interface_Regions ), 1);
        for idx = 1:length( obj.junction.Interface_Regions )
            FD_lead_p{idx} = obj.Fermi(obj.Evec - bias_lead(idx));
            FD_lead_m{idx} = 1 - FD_lead_p{idx};
        end

        % preallocating array
        SelfEnergy = zeros( obj.StatesNum, obj.StatesNum);
       
        
        % calculate the transition rates (only the first order of the perturbation series)
        % from jdx to idx
        % The first order matrix elements are calculated by the diagrammatic rules in PRB 50, 18439
        deltaEvec = [0, diff(obj.Evec)]; 
        for idx = 1:obj.StatesNum             
            
            for jdx = 1:obj.StatesNum
                if idx-1==jdx        
                    % calculating Fermi-Dirac distribution functions in the QD
                    mu_N = obj.mu_QD(jdx);             
                    FD_QD_p = obj.Fermi(obj.Evec-mu_N-obj.Ec);
                    FD_QD_m = 1-FD_QD_p;      
                    
                    % calculating the matrix elements
                    % summation for the leads
                    for kdx = 1:length( obj.junction.Interface_Regions )
                        tmp = 2*1i*obj.TunnelingRates{kdx}(idx,:).*FD_lead_p{kdx}.*FD_QD_m;
                        tmp = obj.SimphsonInt( tmp.*deltaEvec );
                        SelfEnergy(jdx, idx) = SelfEnergy(jdx, idx) + tmp;
                    end
                elseif idx+1 == jdx
                    % calculating Fermi-Dirac distribution functions in the QD
                    mu_N = obj.mu_QD(jdx);
                    FD_QD_p = obj.Fermi(obj.Evec-mu_N);
                    FD_QD_m = 1-FD_QD_p;       
                    
                    % calculating the matrix elements
                    % summation for the leads
                    for kdx = 1:length( obj.junction.Interface_Regions )
                        tmp = 2*1i*obj.TunnelingRates{kdx}(jdx,:).*FD_lead_m{kdx}.*FD_QD_p;
                        tmp = obj.SimphsonInt( tmp.*deltaEvec );
                        SelfEnergy(jdx, idx) = SelfEnergy(jdx, idx) + tmp;
                    end                   
                elseif idx == jdx

                    
                    % calculating the matrix elements
                    % summation for the leads
                    for kdx = 1:length( obj.junction.Interface_Regions )
                        % calculating Fermi-Dirac distribution functions in the QD
                        mu_N = obj.mu_QD(jdx);
                        FD_QD_p = obj.Fermi(obj.Evec-mu_N);
                        FD_QD_m = 1-FD_QD_p;          
                        tmp = -2*1i*obj.TunnelingRates{kdx}(jdx,:).*FD_lead_m{kdx}.*FD_QD_p;
                        
                        if idx < obj.StatesNum
                              mu_N = obj.mu_QD(jdx);
                              FD_QD_p = obj.Fermi(obj.Evec-mu_N-obj.Ec);
                              FD_QD_m = 1-FD_QD_p;
                            tmp = tmp - 2*1i*obj.TunnelingRates{kdx}(idx+1,:).*FD_lead_p{kdx}.*FD_QD_m;
                        end
                        tmp = obj.SimphsonInt( tmp.*deltaEvec );
                        SelfEnergy(jdx, idx) = SelfEnergy(jdx, idx) + tmp;
                    end                      
                else
                   continue
                end               
            end
        end
        
        obj.SelfEnergy = SelfEnergy;
        
        
    end
    
%% CalcTunnelingRates
%> @brief Determine the tunneling rates for all energy values in the energy array Evec for a given gateVoltage
%> @param gateVoltage The back gate voltage in units of the energy unit used in the Hamiltonian parameters.
    function CalcTunnelingRates( obj, gateVoltage )
        
        % create energy array if it is empty
        if isempty( obj.Evec )
            obj.CreateEnergyArray( 'tuned_contact', 'symmetric' );
        end        

        lead_num = length( obj.junction.Interface_Regions );
        
        % preallocating arrays for the energy resolved tunneling rates
        TunnelingRates = zeros( length(obj.Evec), obj.StatesNum, lead_num );    
                
        
    
        % creating handles and temporary variables
        Evec = obj.Evec;
        junction_loc = obj.junction;
        hCalcTunnelingRate = @obj.CalcTunnelingRate;
        
        % starting parallel pool
        poolmanager = Parallel( obj.junction.FL_handles.Read('Opt') );
        poolmanager.openPool(); 
        
       
        % calculating the tunneling rates
        %for idx = 1:length(Evec)
        parfor (idx = 1:length(Evec), poolmanager.NumWorkers) 
            
            Energy = Evec(idx);
            TunnelingRates_tmp = feval( hCalcTunnelingRate, Energy, gateVoltage, 'junction', junction_loc );
            
            % contact resolved tunneling rates
            for iidx = 1:lead_num
                TunnelingRates(idx, :, iidx) = TunnelingRates_tmp(:,iidx);
            end                    
            
        end
        
        %> closing the parallel pool
        poolmanager.closePool();        
        
        obj.TunnelingRates = cell( lead_num, 1);
        for iidx = 1:length( obj.junction.Interface_Regions )
            obj.TunnelingRates{iidx} = transpose(TunnelingRates(:,:,iidx));
        end 
        

%figure
%hold on
%for kdx = 1:obj.StatesNum
%    plot( obj.Evec, real(obj.TunnelingRates{1}(kdx,:)))
%end
        
    end
    
%% CalcTunnelingRate
%> @brief Determine the tunneling rate for a given Energy and gateVoltage 
%> @param Energy The current energy value in units of the energy unit used in the Hamiltonian parameters.   
%> @param gateVoltage The back gate voltage in units of the energy unit used in the Hamiltonian parameters.   
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'junction' An instance of class #NTerminal (or its derived class) representing the junction.
    function TunnelingRates = CalcTunnelingRate( obj, Energy, gateVoltage, varargin )
        
        p = inputParser;
        p.addParameter('junction', obj.junction); %for parallel computations
        p.parse(varargin{:});       
        junction_loc     = p.Results.junction;     
        
        % setting the energy
        obj.SetEnergy( Energy, 'junction', junction_loc ) 
        
        recalculateSurface = 1:length( junction_loc.Interface_Regions );
        
        % preallocating arrays
        TunnelingRates = zeros( obj.StatesNum, length( junction_loc.Interface_Regions ) );
    
        for idx = 1:obj.StatesNum
                        
            % creating function handle for the charging energy effective one-electron potential
            hChargingPotential = @(coordinates)( ones(length(coordinates.x),1)*obj.Ec*( obj.OccNumbers(idx)-gateVoltage/obj.Ec) );
        
            % calculating the surface Green operator of the QD
            junction_loc.CalcFiniteGreensFunctionFromHamiltonian( 'scatterPotential', {obj.scatterPotential, hChargingPotential}, 'onlyGinv', true)            
            
            % Calculate the total Green Operator
            [G, Ginverz, junction_sites] = junction_loc.CustomDysonFunc( ...
                    'decimate', false, 'constant_channels', false, 'recalculateSurface', recalculateSurface, ...
                    'SelfEnergy', obj.useSelfEnergy);                 

            % leads are recalculated only for idx=1
            recalculateSurface = [];
                                    
            % Calculate the spectral funciton
            Spectral = (G - G')/(2*pi*1i);   %check this relation   
            
            % determine the sites corresponding to the leads
            scatter_sites = junction_sites.Scatter.site_indexes;            
            
            for iidx = 1:length(junction_sites.Leads)
                % determine the sites corresponding to the leads
                lead_sites = junction_sites.Leads{iidx}.site_indexes;
                
                % determine the coupling matrix and its adjungate
                coupling = Ginverz( lead_sites, scatter_sites);
                coupling_adj = Ginverz( scatter_sites, lead_sites);
                
                % calculating the tunneling rate
                TunnelingRates(idx, iidx) = trace( coupling_adj* Spectral(lead_sites, lead_sites) * coupling* Spectral(scatter_sites, scatter_sites)); 
                
                
                
            end
            
        end   
        
    end  
    
%% SetChemicalPotentials
%> @brief Manually sets the chemical potentials corresponding to the given single-electron doping levels.
%> @param doping_levels The user defined single-electron doping levels obatined from the calculated density of states.  
%> @param gateVoltage The back gate voltage in units of the energy unit used in the Hamiltonian parameters.   
    function SetChemicalPotentials( obj, doping_levels, gateVoltage )   
        
        % setting the number of energy levels relevant in the calculations
        obj.StatesNum = length(doping_levels);  
        
        % setting the occupation numbers
        % charge neutral occupation is set to be the closest to the zero energy level.
        [~, charge_neutral_index] = min( abs(doping_levels) );
        obj.OccNumbers = (1:obj.StatesNum) -  charge_neutral_index;    
        
        % calculating the chemical potentials icluding the charging energy
        obj.mu_QD = doping_levels + obj.Ec*(obj.OccNumbers - gateVoltage/obj.Ec);
        
        
        %disp(obj.mu_QD)
        
    end
    
%% Initialize
%> @brief Initializes class attributes.
    function Initialize(obj)
        
        obj.Ec = 0;
        obj.StatesNum = 9;
        obj.Edb = 50;
        obj.bias_max = 10*obj.k_B*obj.T;

        if strcmpi(class(obj), 'QuantumDot')             
            obj.InputParsing( obj.varargin{:} );
        end
        
   
         
    end            
    
  
    
    end % public methods
    
    
    methods (Access=protected)    
        
        
%% 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 absolute temperature in Kelvins.
%> @param 'silent' Set true to suppress output messages.
%> @param 'scatterPotential' A function handle y=f( #coordinates coords ) or y=f( #CreateHamiltonians CreateH, E ) of the potential across the junction.
%> @param 'useSelfEnergy' Logical value. Set true (default) to solve the Dyson equation with the self energies of the leads, or false to use the surface Green operators.
%> @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 'junction' An instance of a class #NTerminal or its subclass describing 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 
%> @param 'Ec' The charging energy in the same units as the Hamiltonians.
    function InputParsing( obj, varargin )
        
        p = inputParser;
        p.addParameter('T', obj.T);
        p.addParameter('silent', 0);
        p.addParameter('scatterPotential', []);
        p.addParameter('junction', []);   
        p.addParameter('gfininvfromHamiltonian', false);   
        p.addParameter('useSelfEnergy', true);      
        p.addParameter('Ec', obj.Ec);  
        p.addParameter('bias_max', obj.bias_max);  % The maximal bias to be used in the calculations
        p.addParameter('Edb', obj.Edb);  % The number of the energy points in the attribute Evec
        
        
        p.parse( varargin{:});
        
        InputParsing@Density( obj, 'T', p.Results.T, ...
                            'junction', p.Results.junction, ...
                            'scatterPotential', p.Results.scatterPotential, ...
                            'gfininvfromHamiltonian', p.Results.gfininvfromHamiltonian, ...
                            'useSelfEnergy', p.Results.useSelfEnergy);
         
        obj.Ec                  = p.Results.Ec;
        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