CalculateTransport_pnp_Specq.m

Go to the documentation of this file.

%    Calculates the conductance for specific transverse momentums - based on EQuUs v5.0.144
%    Copyright (C) 2019 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/.
%
%> @file CalculateTransporSpecq.m
%> @brief Function to calculate the transverse momentum resolved conductance through a pnp junction in the zero temperature limit.
%
%
%> @brief Function to calculate the transverse momentum resolved conductance through a pnp junction in the zero temperature limit.
%> @param qvec The transverse momentum number.
%> @param height The length of the junction
%> @param Energy The energy relative to the Fermi level
%> @param Opt An instance of structure #Opt
%> @param param An instance of structure #param
%> @param outputXML The absolute path to the output XML
%> @param pnpPotentialStrength The strength of the npn transition
%> @param B_loc The strength of the magnetic field
%> @return Returns with the calculated conductivity for the given transverse momentum number.
    function Cq = CalculateTransport_pnp_Specq( qvec, height, Energy, EF, Opt, param, outputXML, pnpPotentialStrength, B_loc )
        Cq = zeros(size(qvec));
        
        % Planck contant
        h = 6.626e-34;
        % The charge of the electron
        qe = 1.602e-19;
        % atomic distance
        rCC = 1.42*1e-10; %In Angstrom
        % flux quantum
        phi0 = h/qe;        
        
        % lattice constant in zigzag edged graphene riibon is sqrt(3)*rCC
        lattice_constant = sqrt(3);
        
        % The dimensionless strength of the local magnetic field
        eta_B = 2*pi/phi0*(rCC)^2*B_loc;
        
        for idx = 1:length(qvec)       
            q = qvec(idx);
            
            % shift the transverse momentum with the value of the vector potential to get teh mechanical momentum
            q = q - 1.5*eta_B*height*lattice_constant;
            
            % creating the Ribbon class representing the twoterminal setup
            cRibbon = Ribbon('width', 2, 'height', height, 'Opt', Opt, 'param', param, ...
                    'q', q, 'EF', EF, 'filenameOut', outputXML );
                
            % calculate the conductance for a given transverse momentum
            Cq(idx) = Transport_q( Energy, B_loc, cRibbon );
        end
        
        
    
%% Transport_q
%> @brief Calculates the transvere momemntum resolved conductance for a given energy
%> @param Energy The energy value.
%> @param B The magnetic field
%> @param selfEnergy Logical value. Set true to use the self energy in the calculations or false to use the surface Green operator.
    function Conductance = Transport_q( Energy, B_loc, cRibbon_loc )
        
        % creating funcfion handles for the magnetic vector potentials
     CreateHandlesForMagneticField( B_loc, cRibbon_loc ) 
        
        % seeting the Energy value in the Ribbon class
     cRibbon_loc.setEnergy( Energy )
        
        % Calculates the surface Green operator of the scattering region
        cRibbon_loc.CalcFiniteGreensFunctionFromHamiltonian('PotInScatter', @pnpPotential); 
        
         % creating function handle for the Dyson Eq.
        Dysonfunc = @()cRibbon_loc.CustomDysonFunc( 'constant_channels', false, 'SelfEnergy', true );
        
        % Evaluate the Dyson Eq.
        cRibbon_loc.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
        
        % Calculating the Scattering matrix
        cRibbon_loc.FL_handles.SmatrixCalc();
        
        % Calculate the conductance
        try
          conductance = cRibbon_loc.FL_handles.Conduktance();
        catch
            Conductance = NaN;
            return
        end
        
          Conductance = conductance(1,2);
        
          %disp( ['E = ', num2str(Energy), ' conductance = ', num2str(Conductance)])   
    
    end


%% CreateHandlesForMagneticField
%> @brief Creates and set function handles of the magnetic vector potentials in the Ribbon class
%> @param B The magnetic field
    function CreateHandlesForMagneticField( B_loc, cRibbon_loc )        
        hLandaux = createVectorPotential( B_loc );        
        cRibbon_loc.setHandlesForMagneticField('scatter', hLandaux, 'lead', hLandaux );          
    end
    

%% createVectorPotential
%> @brief Creates the function handle of the magnetic vector potential
%> @param B The magnetic field
%> @return Returns with the function handle. 
    function hLandaux  = createVectorPotential( B_loc )
        
        % lattice constant in zigzag edged graphene riibon is sqrt(3)*rCC
        lattice_constant = sqrt(3);
    
          % calculating the diemnsionless strength of the magnetic field
        eta_B = 2*pi/phi0*(rCC)^2*B_loc;
        % constant vector potential in the systems: stabilizes the numerical computation
        Aconst = 0;         
        
                
        % creting the funciton handles of the vector potentials
        hLandaux = @(x,y)(Landaux(x,y, eta_B, Aconst, height, lattice_constant));
        
    end

%% pnpPotential
%> @brief Potential in the scattering region (transversally translational invariant)
%> @param An instance of structure #coordinates containing the position of the sites
%> @return An array of the on-site potential.
    function ret = pnpPotential( coordinates )
        if isempty(pnpPotentialStrength)
            ret = 0;
            return;
        end
        
        lattice_const = norm(coordinates.a);
        
        lambda = 60*lattice_const; %transition length of the pn potential in units of rCC
        y1         = 0.20*height*lattice_const;
        y2         = height*lattice_const - y1;
        
        y = coordinates.y;        
        
        
         ret = pnpPotentialStrength*( tanh((y-y1)/lambda) - tanh((y-y2)/lambda) )/2;
         
    end
        
        
    end