Graphene_Conductivity.m

Go to the documentation of this file.

%    Energy dependent conductivity through a wide (4 micron) graphene ribbon calculated by the adatpiveQ interface
%    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 Graphene_Conductivity.m
%> @brief Example to calculate the conductivity through a wide (4 micron) graphene junction.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @Available
%> EQuUs v4.8 or later
%> @expectedresult
%> @image html Graphene_Conductivity.jpg
%> @image latex Graphene_Conductivity.jpg
%> The red line represents the theoretical \f$ 2/\pi\;\sigma_0 \f$ limit of the minimal conductivity with \f$ \sigma_0 = e^2/\hbar \f$.
%> @}
%
%> @brief Example to calculate the conductivity through a wide (4 micron) graphene junction
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function Graphene_Conductivity( 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 data
    outfilename = [fncname,'_',num2str( filenum )];
    % the output folder
    outputdir   = [];

    % length of the junction
    height = 2400;
    % width of the junction 
    total_width = 19200;
    % width of the unit cell
    width_uc = 2; 
    % The aspect ratio of the junction
    aspect_ratio = (total_width*3/2)/((height)*sqrt(3));

    % The Fermi energy in the calculations
    EF = 0; %In EV
    % 1D energy array contaning energy values in the transport calculations
    EF_vec = [1e-5:1e-3:0.04];
    % Energy ... used in transport calculations with full Hamiltonian
    DeltaEF_index = 5;
    % 1D energy array contaning energy values in the transport calculations with full Hamiltonian
    EF_vec_fromH = EF_vec(1:DeltaEF_index:end);
    % Energy in the transport calculations relative to the Fermi energy
    Energy = 0;

    % Conductivity calculated by the fast way
    Conductivity = zeros(size(EF_vec));

    % conductivity calculated from the full Hamiltonian
    Conductivity_fromH = zeros(size(EF_vec_fromH));




    % creating output directories
    setOutputDir();    
    
    % filename containing the output XML
    outputXML = [outputdir, filesep, outfilename, '.xml'];
    
    % Calculating the conductivity
    disp(['Performing calculation for junction with  aspect ratio W/L = ', num2str(aspect_ratio, '%10.2f')]);
    calcConductivity();
    
    % export the calculetd data
    save( [outputdir,'/',outfilename,'.mat']);  
    
    % nested functions

    %% calcConductivity
    %> @brief Calculate the conductivity for the selected energy points
    function calcConductivity
        
        % for graphene infinite mass boundary condition is the zigzag orientated ibbon, q in units of 1/(3r_{cc})
        W = (total_width/2*2 + (total_width/2-1) ); %in units of r_CC
        m = 0:1:total_width/width_uc;
        interpq = transpose( (m+0.5)*pi/(W/3) ); % for graphene infinite mass boundary
        
                
        for idx = 1:length(EF_vec)
        
            disp(' ')
            disp(' ')
            disp(' ----------------------------------')
            disp('Calculating Conductivity for the given Fermi energy')    
            
            EF = EF_vec(idx);
            
            % function handle for the transport calculation using the fast method
            TransportFast = @(qvec)(CalculateTransporSpecq( qvec, width_uc, height, Energy, EF, Opt, param, outputXML, false ) );
            
            % creating an instance of class adaptiveQ for transverse momentum iterations
            adaptive_moments = adaptiveQ( Opt, interpq );
            
            % running the transverse momentum iterations with the function handle for the transport calculations
            adaptive_moments.runAdaptiveIterations( TransportFast );
            
            % collecting the calculated data
            sumCq = sum(adaptive_moments.Read('interpy'))/aspect_ratio;
            Conductivity(idx) = sumCq;
            
            if mod(idx-1,DeltaEF_index) == 0
                disp(' ')
                disp(' ')
                disp(' ----------------------------------')
                disp('Calculating Conductivity for the given Fermi energy from the full Hamiltonian')
                
                % function handle for the transport calculation from the full Hamiltonian
                TransportSlow = @(qvec)(CalculateTransporSpecq( qvec, width_uc, height, Energy, EF, Opt, param, outputXML, true) );
                
                % creating an instance of class adaptiveQ for transverse momentum iterations
                adaptive_moments = adaptiveQ( Opt, interpq );
                
                % running the transverse momentum iterations with the function handle for the transport calculations
                adaptive_moments.runAdaptiveIterations( TransportSlow );
                
                % collecting the calculated data
                sumCq = sum(adaptive_moments.Read('interpy'))/aspect_ratio;
                Conductivity_fromH( floor(idx/DeltaEF_index)+1 ) = sumCq;
            end
            
            disp(' ----------------------------------')
            disp('Plotting conductivity:')
            PlotFunction();        
            
            disp([ num2str(idx/length(EF_vec)*100),' % calculated of the conductance.'])
            save( fullfile(outputdir, [outfilename, '.mat']), 'Conductivity', 'Conductivity_fromH', 'Energy', 'EF_vec', 'height', 'total_width');
        end
        
        
        

    end











%% setOutputDir
%> @brief Sets the output directory
    function setOutputDir()
        resultsdir = fullfile(pwd, 'results');
        mkdir(resultsdir );
        outputdir = resultsdir;        
    end

%% PlotFunction
%> @brief Plots the calculated data
    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;
        
        %Theoretical minimal conductivity
        theoretical_limit = 2/pi; 
        
        % determine points to be plotted
        indexek = logical(Conductivity);
        indexek_fromH = logical(Conductivity_fromH);
        x_lim = [0 max(EF_vec(indexek))*1000];
        
        if norm(x_lim) == 0
            x_lim = [0 1];
        end
        
        % creating axes of the plot   
        axes_cond = axes('Parent',figure1,... 'Position', Position,...
                'Visible', 'on',...
                'FontSize', fontsize,... 
                'xlim', x_lim,...  
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        legend_labels = [];
        
        % plot the data
        if ~isempty(Conductivity(indexek))
            numerics = plot(EF_vec(indexek)*1000, Conductivity(indexek), 'Linewidth', 2, 'color', [0 0 0], 'Parent', axes_cond);  
            legend_labels{end+1} = 'PRB 90, 125428 (2014)';
        end
        if ~isempty(Conductivity_fromH(indexek_fromH))
            numerics2 = plot(EF_vec_fromH(indexek_fromH)*1000, Conductivity_fromH(indexek_fromH), 'LineStyle', 'none', 'Marker', 'v', 'MarkerSize',6, 'color', [0 0 1], 'Parent', axes_cond);   
            legend_labels{end+1} =  'EQuUs MKL';
        end
        
         plot( [0 max(EF_vec(indexek))*1000], ones(1,2)*theoretical_limit, 'r',  'Parent', axes_cond);
        
               
        % Create xlabel
        xlabel('E_F [meV]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_cond);

        % Create ylabel
        ylabel('\sigma/\sigma_0','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_cond);
        
        if length( legend_labels) == 2
            fig_legend = legend([numerics, numerics2], legend_labels);
            set(fig_legend, 'FontSize', fontsize, 'FontName','Times New Roman', 'Box', 'off', 'Location', 'NorthEast')
        end
        
        % 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, [fncname, '_', num2str(filenum),'.eps']))
        print('-dpdf', fullfile(outputdir, [fncname, '_', num2str(filenum),'.pdf']))
        close(figure1)
        
    end


    % end nested functions

end