magnetic_barrier.m

Go to the documentation of this file.

%   Transport through magnetic barrier
%    Copyright (C) 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/.
%
%> @addtogroup Examples
%> @{
%> @file magnetic_barrier.m
%> @brief Example to calculate the conductivity through a magnetic barrier.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @Available
%> EQuUs v4.8 or later
%> @expectedresult
%> @image html magnetic_barrier.jpg
%> @image latex magnetic_barrier.jpg
%> @}
%
%> @brief Example to calculate the conductivity through a magnetic barrier.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function magnetic_barrier( filenum )

if ~exist('filenum', 'var')
    filenum = 1;
end

filename = mfilename('fullpath');
[directory, fncname] = fileparts( filename );

    % filename containing the input XML    
    inputXML = 'Basic_Input_zigzag_leads.xml';
    % Parsing the input file and creating data structures
    [Opt, param] = parseInput( fullfile( directory, inputXML ) ); 
    
    % filename containing the output XML
    outfilename = [fncname, '_',num2str( filenum )];
    % the output folder
    outputdir   = [];    
    
    % width of the junction (number of non-singular sites in the cross section)
    width = 160;
    % length of the junction (number of unit cells)
    height = 300;    

    % creating 1D energy array (in units of eV)
    % minimum of the energy array
    Emin = 0.001;
    % maximum of the energy array
    Emax = 0.6;
    %> number of energy points
    Enum = 200;
    % the energy array
    Evec = Emin:(Emax-Emin)/(Enum+1):Emax;

    % the strength of the magnetic field in Tesla
    B = 40e-0;
    % the cyclotron frequency
    homegac = [];
    % class representing the two terminal ribbon structure
    cRibbon = [];
    % Planck contant
    h = 6.626e-34;
    % The charge of the electron
    qe = 1.602e-19;
    % atomic distance
    rCC = 1.42*1e-10; %In Angstrom

    % preallocating the array for the calculated conductance values
    Conductance = zeros(1,length(Evec)); 
    % preallocating the array for the unitarity error
    DeltaC      = zeros(1,length(Evec));

    % creating output directories
    setOutputDir();



    % Calculate the transport through the magnetic barrier
    CalculateTransport();
    
    % plot the calculated results
    PlotFunction();




%% CalculateTransport
%> @brief Calculates the conductance
    function CalculateTransport()
        
        % find indexes needed to be calculated
        indexes2calc = find( Conductance == 0 );
        
        
        % perform the transport calculations
        for idx = 1:length(indexes2calc)  
            
            % determine the current energy value
            Energy = Evec(indexes2calc(idx));    
            
            % creating the Ribbon class representing the twoterminal setup
            cRibbon = Ribbon('Opt', Opt, 'param', param, 'width', width, 'height', height );
            
            % calculating dimensionless magnetic field, and the cyclotron frequency
            phi0 = h/qe; % magnetic flux quantum
            eta_B = 2*pi/phi0*(rCC)^2*B;
            if B>1e-8
                homegac = 2.97*sqrt(9/2*eta_B);
            end

            % Plot the calculated results
            disp(' ----------------------------------')
            disp('Plotting conductivity:')
            PlotFunction();        
        
        
            % calculating the current conductace
            disp(' ')
            disp(' ')
            disp(' ----------------------------------')
            disp('Calculating Conductance for the given magentic field')  

            
            tic
            [Conductance_tmp,DeltaC_tmp] = Transport( Energy, B );
            toc
            
            % storing the calculated data
            Conductance(indexes2calc(idx)) = Conductance_tmp;
            DeltaC(indexes2calc(idx)) = DeltaC_tmp;
            
            
            % exporting the calculated data
            save( fullfile(outputdir, [outfilename, '.mat']));
            disp([ num2str(indexes2calc(idx)/length(Evec)*100),' % calculated of the conductance.'])
            disp( ' ' )
            
            
            
        end
                
    end

%% Transport
%> @brief Calculates the conductance at a given energy value.
%> @param Energy The energy value.
%> @param B The magnetic field
%> @return Returns with teh calculated conductance and the unitarity error.
    function [Conductance,DeltaC] = Transport( Energy, B )
    
        % creating funcfion handles for the magnetic vector potentials
     CreateHandlesForMagneticField( B )      
        
        % seeting the Energy value in the Ribbon class
     cRibbon.setEnergy( Energy )
        % Calculates the surface Green operator of the scattering region and perform a gauge transformation
        cRibbon.CalcFiniteGreensFunction( 'gauge_trans', true );  
        
        % creating function handle for the Dyson Eq.
        Dysonfunc = @()cRibbon.CustomDysonFunc( 'constant_channels', true, 'SelfEnergy', false );
        
        % Evaluate the Dyson Eq.
        cRibbon.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
        
        % Calculating the Scattering matrix
        S = cRibbon.FL_handles.SmatrixCalc();
        
        % Calculate the unitarity error
        DeltaC = norm(S*S'-eye(size(S)));
        
        if DeltaC >= 1e-3 || isnan(DeltaC)
            disp( ['error of the unitarity of S-matrix: ',num2str(DeltaC)] )
        end
        
        
        % Calculate the conductance
        try
          Conductance = cRibbon.FL_handles.Conduktance();
            Conductance = Conductance(1,2);
        catch
            Conductance = NaN;
            DeltaC = NaN;
            
            return
        end
        
          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 ) 
        [hLead, hScatter, gauge_field ] = createVectorPotential( B );        
        cRibbon.setHandlesForMagneticField('scatter', hScatter, 'lead', hLead, 'gauge_field', gauge_field );
    end
    

%% createVectorPotential
%> @brief Creates the function handles of the magnetic vector potentials and the gauge field
%> @param B The magnetic field
%> @return [1] Function handle of the vector potential in the leads .
%> @return [2] Function handle of the vector potential in the scattering region.
%> @return [3] Function handle of the scalar potential that transforms vector potential in the scattering center into a vector potential in the leads.
    function [hLead, hScatter, hgauge_field ] = createVectorPotential( B )
        
          phi0 = h/qe;
          
        % calculating the diemnsionless strength of the magnetic field
        eta_B = 2*pi/phi0*(rCC)^2*B;
        % constant vector potential in the systems: stabilizes the numerical computation
        Aconst = 0.1;  
        
        % lattice constant in zigzag edged graphene riibon is sqrt(3)*rCC
        lattice_constant = sqrt(3);
                
        % creting the funciton handles of the vector potentials
        hLead = @(x,y)(Landaux(x,y, eta_B, Aconst, height, lattice_constant));
        hScatter = @(x,y)(Landauy(x,y, eta_B));
        hgauge_field = @(x,y)(gauge_field(x,y, eta_B, Aconst, height, lattice_constant));        
        
    end


%% setOutputDir
%> @brief Sets output directory.
    function setOutputDir()
        resultsdir = 'results';
        mkdir(resultsdir );
        outputdir = [ resultsdir ];        
    end



%% PlotFunction
%> @brief Creates the plot
    function PlotFunction( )
        

        % creating figure in units of pixels
        figure1 = figure( 'Units', 'Pixels', 'Visible', 'off');
        
        % font size on the figure will be 16 points
        fontsize = 16;
                       
        % determine points to be plotted
        indexek = logical( (abs(DeltaC) <= 0.5) & (Evec~=0) & ~isnan(DeltaC) );
        
        % setting the limits of the x axis
        x_lim = [min(Evec(indexek)) max(Evec(indexek))];        
        
        if norm(x_lim) == 0
            x_lim = [0 1];
        end
        
        % creating axes of the plot        
        axes_cond = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,... 
                'xlim', x_lim,...            
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        % plot the data
        plot(Evec(indexek), Conductance(indexek), 'Linewidth', 2, 'color', [0 0 0], 'Parent', axes_cond);  
        plot(Evec(indexek), DeltaC(indexek), 'Linewidth', 1, 'color', [0 1 0], 'Parent', axes_cond);
        
        % setting the cyclotron frequencies in the plot
        if ~isempty(homegac)
            nmax = round( (Emax/homegac)^2 );
            nvec = 1:nmax;
            for idx = 1:length(nvec)
                plot(sqrt([nvec(idx), nvec(idx)])*homegac, [0 max(Conductance(indexek))], 'Linewidth', 0.5, 'color', [0 0 1], 'Parent', axes_cond);
            end
        end
        
        
        % Create xlabel
        xlabel('E [ev]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_cond);

        % Create ylabel
        ylabel('\sigma/\sigma_0','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_cond);
        
        % set the legends
        fig_legend = legend(axes_cond, {'conductance', 'unitarity error'});%, 'FontSize', fontsize, 'FontName','Times New Roman')
        set(fig_legend, 'FontSize', fontsize, 'FontName','Times New Roman', 'Box', 'off', 'Location', 'NorthEast')
        
        % setting the position and margins of the plot, removing white
        % spaces for release dates greater than 2015
        ver = version('-release');
        if str2num(ver(1:4)) >= 2016
            figure1.PaperPositionMode = 'auto';
            fig_pos = figure1.PaperPosition;
            figure1.PaperSize = [fig_pos(3) fig_pos(4)]; 
        
            set(axes_cond, 'Position', get(axes_cond, 'OuterPosition') - get(axes_cond, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1]);        
            Position_figure = get(axes_cond, 'OuterPosition');
            set(figure1, 'Position', Position_figure);
        end

        % export the figures
        print('-depsc2', fullfile(outputdir,[outfilename, '.eps']))
        print('-dpdf', fullfile(outputdir,[outfilename, '.pdf']))
        close(figure1);
        
        
    end




end