magnetic_ribbon.m

Go to the documentation of this file.

%   Transport through graphene ribbon in magnetic field -  Based on EQuUs v4.7
%    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 Examples
%> @{
%> @file magnetic_ribbon.m
%> @brief Example to calculate the conductivity through a ribbon in a magnetic field.
%> @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_ribbon.jpg
%> @image latex magnetic_ribbon.jpg
%> @}
%
%
%> @brief Example to calculate the conductivity through a ribbon in a magnetic field.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function magnetic_ribbon( 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   = [];      

    %> length of the junction (number of unit cells)
    height = 300;
    % width of the junction 
    width = 200;

    % 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; %in Tesla
    % the cyclotron frequency
    homegac = [];

    % 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 errors
    DeltaC      = zeros(1,length(Evec));
    % preallocating the array for the difference between the calculations using the self energies and the surface Green function
    Difference  = zeros(1,length(Evec));
    % class representing the two terminal ribbon structure
    cRibbon = [];

    %> creating output directories
    setOutputDir();

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




%% CalculateTranspor
%> @brief Calculates the conductance values
    function CalculateTranspor()
        
        for idx = 1:length(Evec)  
            % determine the current energy value
            Energy = Evec(idx);  
            
            %> creating the Ribbon class representing the twoterminal setup
            cRibbon = Ribbon('width', width, 'height', height, 'Opt', Opt, 'param', param, 'filenameOut', fullfile( outputdir, [outfilename, '.xml']) );
            
            % 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, false );
            [Conductance_tmp_selfEnergy] = Transport( Energy, B, true );
            toc
            
            % storing the calculated data
            Conductance(idx) = Conductance_tmp;
            DeltaC(idx) = DeltaC_tmp;
            Difference(idx) = Conductance_tmp - Conductance_tmp_selfEnergy;
            
            
            % exporting the calculated data
            save( fullfile(outputdir,[outfilename, '.mat']) );
            disp([ num2str(idx/length(Evec)*100),' % calculated of the conductance.'])
            disp( ' ' )
            
        end
                
    end
    
%% Transport
%> @brief Calculates the 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,DeltaC] = Transport( Energy, B, selfEnergy )
        
        % 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
        cRibbon.CalcFiniteGreensFunction(); 
        
         % creating function handle for the Dyson Eq.
        Dysonfunc = @()cRibbon.CustomDysonFunc( 'constant_channels', false, 'SelfEnergy', selfEnergy );
        
        % Evaluate the Dyson Eq.
        cRibbon.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
        
        % Calculating the Scattering matrix
        [S,open_channels] = cRibbon.FL_handles.SmatrixCalc();
        
        % Calculate the unitarity error
        DeltaC = norm(S*S'-eye(sum(open_channels)));
        
        if DeltaC >= 1e-3 || isnan(DeltaC)
            disp( ['error of the unitarity of S-matrix: ',num2str(DeltaC)] )
        end
        
        if open_channels(1) ~= open_channels(2)
           disp( 'openchannels do not match' )
        end
        
        % Calculate the conductance
        try
          conductance = cRibbon.FL_handles.Conduktance();
        catch
            Conductance = NaN;
            DeltaC = 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 )        
        hLandauy = createVectorPotential( B );        
        cRibbon.setHandlesForMagneticField('scatter', hLandauy, 'lead', hLandauy );          
    end
    

%% createVectorPotential
%> @brief Creates the function handle of the magnetic vector potential
%> @param B The magnetic field
%> @return Returns with the function handle. 
    function hLandauy  = createVectorPotential( B )
    
          phi0 = h/qe;        
        eta_B = 2*pi/phi0*(rCC)^2*B;
        
        % creting the funciton handles of the vector potential
        hLandauy = @(x,y)([zeros(size(x,1),1),eta_B*x]);
        
    end



%% setOutputDir
%> @brief Sets output directory.
    function setOutputDir()
        resultsdir = 'results';
        mkdir(resultsdir );
        outputdir = fullfile( directory, 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
        numerics = plot(Evec(indexek), Conductance(indexek), 'Linewidth', 2, 'color', [0 0 0], 'Parent', axes_cond);  
        try
            plot(Evec(indexek), DeltaC(indexek), 'Linewidth', 1, 'color', [0 1 0], 'Parent', axes_cond);
            plot(Evec(indexek), Difference(indexek), 'Linewidth', 1, 'color', [1 0 0], 'Parent', axes_cond);
        catch
            disp('No DeltaC data present')
        end
        
        % 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', '"Surface Green Func - Self energy"'});
        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
            fig = gcf;
            fig.PaperPositionMode = 'auto';
            fig_pos = fig.PaperPosition;
            fig.PaperSize = [fig_pos(3) fig_pos(4)]; 
        
            set(gca, 'Position', get(gca, 'OuterPosition') - get(gca, '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