SpectralDensity.m

Go to the documentation of this file.

%   Spectral density function in superconductor-graphene-superconductor junctions.
%    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 SpectralDensity.m
%> @brief Example to calculate the spectral density function in a superconductor-graphene-superconductor junction.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @Available
%> EQuUs v4.7 or later
%> @expectedresult
%> @image html SpectralDensity.jpg
%> @image latex SpectralDensity.jpg
%> @}
%
%> @brief Example to calculate the spectral density function in a superconductor-graphene-superconductor junction.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%%
function SpectralDensity( 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   = [];     

     % The width of the unit cell
    width = 2;
     % The length of the graphene strip
    height = 1016;
    
    % setting the output directory
     setOutputDir()
    
     % outfilename to export the running parameters
     outputXML = fullfile(outputdir, [outfilename, '.xml']);    

    % carbon-carbon atom distance
    rCC = 1.42e-10;
    % the phase difference between the superconducting contacts
    deltaPhi = 1.4043;
    % The Fermi energy in the scattering region
    EF = 0.08; %In EV
    % imaginary part of the Fermi energy
    eta = 1e-6;
    % parameter for the doping profile
    surfacePotential = 0.1; %in eV
    
    % setting the phase of the superconducting pair potentials
    param.Leads{1}.pair_potential = abs(param.Leads{1}.pair_potential);
    param.Leads{2}.pair_potential = abs(param.Leads{2}.pair_potential)*exp(1i*deltaPhi);
    
    
    % creating array containing the pair potentials in the leads
    Delta = [param.Leads{1}.pair_potential, param.Leads{2}.pair_potential ];
    
    % setting the energy and tranverse momentum grids for the calculations
    [Evec, qvec] = setVectors();

    % An instance of class Ribbon describing the junction
    cRibbon    = [];
    
    % The calculated density of states
    density_of_states = zeros( length(Evec), length(qvec));

    % calculate the density of states
    CalculateSpectralDensity();




%% CalculateSpectralDensity
%> @brief Calculate the density of states
    function CalculateSpectralDensity()
        
        % determine indexes to be calculated
        indexes = logical( density_of_states(1,:) );
        index2calc = 1:length(indexes);
        index2calc = index2calc( ~indexes );
        
        
        % Do the calculations
        for idx = 1:length(index2calc)
            q = qvec(index2calc(idx));
            
            disp(' ')
            disp(' ')
            disp(' ----------------------------------')
            disp(['Calculating Greens function for the transverse momentum q=', num2str(q)]) 
                
            % temporary array to store the calculated DOS
            Rho = zeros( length(Evec), 1);
            
            % Creating an instance of Ribbon class with the current transverse momentum q
            cRibbon = Ribbon('width', width, 'height', height, 'E', 0, 'EF', EF, ...
                'Opt', Opt, 'param', param, 'filenameOut', outputXML, 'q', q);
                
            for jdx=1:length(Evec)
                Energy = Evec(jdx);
                
                % setting the current energy
                cRibbon.setEnergy( Energy+1i*eta );
                
                % Calculates the DOS for the given energy and transverse momentum quantum number
                [Atot, G] = cRibbon.CalcSpectralFunction( Energy+1i*eta, 'constant_channels',false, ...
                'gfininvfromHamiltonian', true, 'PotInScatter', @SurfacePot, ...
                'decimateDyson', false);
                            
                Rho(jdx) = Atot;                   
                
            end
            
            density_of_states(:,index2calc(idx)) = Rho;
            
            % calculated data are plotted at each step
            PlotFunction( )
            
            % exporting the calculated data
            save( [outputdir,'/',outfilename, '.mat'], 'qvec', 'Evec', 'width', 'height', 'EF', 'density_of_states' );
            
            % displaying the status of the calculations
            disp([ num2str(index2calc(idx)/length(qvec)*100),' % calculated of the density of states.'])
        end
        
        
                
    end


%% SetVectors
%> @brief Sets the energy and transverse quantum number grid points
%> @return [1] 1D array of the energy points.
%> @return [2] 1D array of the transverse momentum points.
    function [Evec, qvec] = setVectors()
        Evec = (0:max(abs(Delta))/70:1.1*max(abs(Delta)));
        vargamma = param.scatter.vargamma;
        qmax = (abs(EF)+max(abs(Evec)))/(vargamma/2)*1.7;
        qvec = (-qmax:qmax/70:qmax); % in units of (3r_CC)^-1

    end

%% ScatterPotetial
%> @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 = SurfacePot( coordinates )
        if isempty(surfacePotential)
            ret = 0;
            return;
        end
        
        lattice_const = norm(coordinates.a);
        
        lambda = 15*lattice_const; %transition length of the pn potential in units of rCC
        y1         = 0.15*height*lattice_const;
        y2         = height*lattice_const - y1;
        
        y = coordinates.y;        
        
        
         ret = surfacePotential*( tanh((y-y1)/lambda) - tanh((y-y2)/lambda) )/2;
    end

%% setOutputDir
%> @brief sets the output directory
    function setOutputDir()
        resultsdir = 'results';
        mkdir(resultsdir );
        outputdir = resultsdir;        
    end

%% PlotFunction
%> @brief plot the calculated data
function PlotFunction( )
 
% ********************** plot the DOS ***********************

        % determine points to be plotted
        indexes = logical( density_of_states(1,:));
        
        if length(qvec(indexes)) < 2
            return
        end
        
        % creating figure in units of pixels
        figure1 = figure( 'Units', 'Pixels', 'Visible', 'off', 'Colormap', ...
    [0.0416666679084301 0 0;0.0833333358168602 0 0;0.125 0 0;0.16666667163372 0 0;0.20833332836628 0 0;0.25 0 0;0.291666656732559 0 0;0.333333343267441 0 0;0.375 0 0;0.416666656732559 0 0;0.458333343267441 0 0;0.5 0 0;0.541666686534882 0 0;0.583333313465118 0 0;0.625 0 0;0.666666686534882 0 0;0.708333313465118 0 0;0.75 0 0;0.791666686534882 0 0;0.833333313465118 0 0;0.875 0 0;0.916666686534882 0 0;0.958333313465118 0 0;1 0 0;1 0.0416666679084301 0;1 0.0833333358168602 0;1 0.125 0;1 0.16666667163372 0;1 0.20833332836628 0;1 0.25 0;1 0.291666656732559 0;1 0.333333343267441 0;1 0.375 0;1 0.416666656732559 0;1 0.458333343267441 0;1 0.5 0;1 0.541666686534882 0;1 0.583333313465118 0;1 0.625 0;1 0.666666686534882 0;1 0.708333313465118 0;1 0.75 0;1 0.791666686534882 0;1 0.833333313465118 0;1 0.875 0;1 0.916666686534882 0;1 0.958333313465118 0;1 1 0;1 1 0.0625;1 1 0.125;1 1 0.1875;1 1 0.25;1 1 0.3125;1 1 0.375;1 1 0.4375;1 1 0.5;1 1 0.5625;1 1 0.625;1 1 0.6875;1 1 0.75;1 1 0.8125;1 1 0.875;1 1 0.9375;1 1 1]);
        
        % font size on the figure will be 16 points
        fontsize = 12;        
        
        % define the colorbar limits
        colbarlimits = [min(min(log(density_of_states))) max(max(log(density_of_states)))];
        
        % define the axis limits
        x_lim = [min(qvec) max(qvec)];
        y_lim = [min(Evec) max(Evec)]/min(abs(Delta));
        
        
        axes_DOS = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,... 
                'xlim', x_lim,...           
                'ylim', y_lim,...      
                'CLim', colbarlimits, ...
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        % plot the data
        [X,Y] = meshgrid( qvec(indexes), Evec );
        levelnum = 50;
        density_of_states2plot = log(density_of_states(:,indexes));
        db=1;
        contour(X(1:db:end,1:db:end), Y(1:db:end,1:db:end)/min(abs(Delta)), density_of_states2plot(1:db:end,1:db:end), levelnum ,'LineStyle','none','Fill','on',...   'LevelList',[0.2391 0.2866 0.3342 0.3817 0.4293 0.4769 0.5244 0.572 0.6195 0.6671 0.7147 0.7622 0.8098 0.8573 0.9049 0.9524],...
           'Parent', axes_DOS);
        
        % Create xlabel
        xlabel('q [1/(3r_{CC})]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_DOS);
        
        % Create ylabel
        ylabel('E [\Delta]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_DOS);



     

% ********************** plot potential ************************
        coordinates = structures('coordinates');
        coordinates.y = 0:0.01:height*sqrt(3);
        coordinates.a = [0, sqrt(3)];
        
        % determine the potential across the svócattering region
          pot = SurfacePot( coordinates );
                 
        % position dependent local doping
        lead_shift = param.Leads{1}.epsilon;
          mu = EF - lead_shift - pot;
        
          x_lim = [ min(coordinates.y-fix(height/4)*coordinates.a(2)) max(coordinates.y+fix(height/4)*coordinates.a(2))]*rCC*1e10;
          y_lim = [ min(mu) max(mu)];
        
        axes_pot = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,... 
                'xlim', x_lim,...   
                   'ylim', y_lim,...    
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        
        % plot the data
          plot( coordinates.y*rCC*1e10, mu, 'LineWidth', 2, 'parent', axes_pot);
        plot( [min(coordinates.y)*rCC*1e10, min(coordinates.y)*rCC*1e10], y_lim, 'LineStyle', '--', 'color', [0 0 0], 'parent', axes_pot )
        plot( [max(coordinates.y)*rCC*1e10, max(coordinates.y)*rCC*1e10], y_lim, 'LineStyle', '--', 'color', [0 0 0], 'parent', axes_pot )
                   

        % setting the xlabel and ylabel
          xlabel('x [A]', 'FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_pot);
          ylabel('doping level [ev]' ,'FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_pot);
        
        % det the position of the figure
        figure_pos = get( figure1, 'Position' );
        figure_pos(4) = figure_pos(4)/2;
        set( figure1, 'Position', figure_pos);
        
        
        %set the position of the axes_DOS
        OuterPosition = get(axes_DOS, 'OuterPosition');
        OuterPosition(1) = 0;
        OuterPosition(2) = 0;
        OuterPosition(3) = figure_pos(3)/2;
        OuterPosition(4) = figure_pos(4);
        set(axes_DOS, 'OuterPosition', OuterPosition);      
        Position_DOS = get(axes_DOS, 'OuterPosition') - get(axes_DOS, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1];
        %set(axes_DOS, 'Position', Position_DOS);              
        
        %set the position of the axes_pot
        OuterPosition = get(axes_pot, 'OuterPosition');
        OuterPosition(1) = figure_pos(3)/2;
        OuterPosition(2) = 0;
        OuterPosition(3) = figure_pos(3)/2;
        OuterPosition(4) = figure_pos(4);
        set(axes_pot, 'OuterPosition', OuterPosition);      
        Position_pot = get(axes_pot, 'OuterPosition') - get(axes_pot, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1];
        %set(axes_pot, 'Position', Position_pot);              
                
        % set the paper position in releases greater than 2016
        ver = version('-release');
        if str2num(ver(1:4)) >= 2016
            % setting the position and margins of the plot, removing white spaces
            figure1.PaperPositionMode = 'auto';
            fig_pos = figure1.PaperPosition;
            figure1.PaperSize = [fig_pos(3) fig_pos(4)/2];  
        end
            
                
        % export the figures
        try
            print('-depsc2', fullfile(outputdir, [outfilename,'.eps']))
            print('-dpdf', fullfile(outputdir,[outfilename, '.pdf']))
        catch errCause
            disp('Failed to export figure');
            disp( errCause.identifier );
            disp( errCause.stack(1) );
        end
        close(figure1)
        
        
    end

   



end