NS_NEGF.m

Go to the documentation of this file.

%%   Energy resolved current density across a biased normal-superconductor interface 
%    Copyright (C) 2018 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 NS_NEGF.m
%> @brief Example to calculate the energy resolved current density across a biased normal-superconductor interface 
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @Available
%> EQuUs v4.9 or later
%> @expectedresult
%> @image html NS_NEGF.jpg
%> @image latex NS_NEGF.jpg
%> @}
%
%> @brief Example to calculate the energy resolved current density across a biased normal-superconductor interface 
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function NS_NEGF( filenum )

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

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

% filename containing the output XML
outfilename = [fncname, '_',num2str( filenum )];
% the output folder
outputdir   = [];   

% filename of the input parametersubuntu download usb installer
filenameIn  = [fncname, '.xml'];

% The width of the junction
width = 40;

% The length of the junction
height = 4;

% number of energy points on the energy contour
Edb = 111;
% the superconducting pair potential        
pair_potential = 0.1;

% loading data structures from the input XML
[Opt, param] = parseInput( fullfile( pwd, filenameIn) );

% set the width and length of the junction
setShape()

% creating the output directory
setOutputDir()

% the energy array
Evec = [];
% the calculated current through the contacts 
currentsurf = [];

% calculate the Josephson current
CalculateNSCurrent()

% Plot the calculated results
PlotFunction( );


%% setWidth
%> @brief setting the width and length of the junction
    function setShape()
        
        for idx = 1:length( param.Leads )
            param.Leads{idx}.M = width;
        end
        
        param.scatter.shape.width = width;
        param.scatter.shape.height = height;
        
    end

%% CalculateJosephson
    function CalculateNSCurrent()
                          
        % turning on the BdG model
        Opt.BdG = 1;
        
        % bias leads
        bias_leads = [0.00 0.15];
        
        % chemical potential in the central device
        EF = 0;
        
        % set new temperature value
        T = 0;  
        
        % setting the superconducting pair potentials in the leads
        param.Leads{1}.pair_potential = pair_potential;
        param.Leads{2}.pair_potential = 0;
        
        param.scatter.pair_potential = pair_potential;
        
        % creating the NTerminal class
        cRibbon_K = Ribbon_Keldysh('Opt', Opt, 'param', param, 'filenameOut', fullfile( outputdir, [outfilename, 'xml']), ...
            'EF', EF, 'bias_leads', bias_leads, 'T', T);           
        
        % creating class to calculate the Josephson effect
        cCurrent_NEGF = Current_NEGF(Opt, 'junction', cRibbon_K, 'T', T, 'gfininvfromHamiltonian', true );
               
        % calculating the Josephson current          
        DeltaPhi_vec = 0;
        [currentvec, currentsurf, Evec] = cCurrent_NEGF.CurrentCalc_discrete( DeltaPhi_vec, 'Emin', EF-0.9*pair_potential, 'Emax', EF+1.45*pair_potential, 'Edb', Edb );
            
                               
    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;
         
        %************** conductance **********************
        
        % creating axes of the plot        
        axes_cond = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,...    
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        % plot the data
        plot(Evec*1000, real(currentsurf(:,2)), 'Linewidth', 2, 'color', [0 0 0], 'Parent', axes_cond);  
        
          
        
        % Create xlabel
        xlabel('$$E\;[meV]$$','FontSize', fontsize, 'Interpreter', 'latex', 'FontName','Times New Roman', 'Parent', axes_cond);

        % Create ylabel
        ylabel('$$I\;[2e^2\Delta/\hbar]$$','FontSize', fontsize, 'Interpreter', 'latex', 'FontName','Times New Roman', 'Parent', axes_cond);
        
        figure_pos = get( figure1, 'Position' );
        
        %set the position of the axis
        OuterPosition = get(axes_cond, 'OuterPosition');
        OuterPosition(1) = 0;
        OuterPosition(2) = 0;
        OuterPosition(3) = figure_pos(3);
        OuterPosition(4) = figure_pos(4);
        set(axes_cond, 'OuterPosition', OuterPosition);      
        Position_cond = get(axes_cond, 'OuterPosition') - get(axes_cond, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1];
        set(axes_cond, 'Position', Position_cond);              
        
        % 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
            % 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)];                  
        end
               

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







%% sets the output directory
    function setOutputDir()
        
        % determine the directory for the results
        resultsdir = fullfile(pwd, 'results');
        % determine the working directory
        workingdir = fullfile(pwd);

        % creating the working directory
        mkdir(workingdir);
        % creating the results directory
        mkdir(resultsdir );
        outputdir = resultsdir;        
        
    end




    end