JosephsonCurrent.m

Go to the documentation of this file.

%    DC Josephson current through a one-dimensional chain
%    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 JosephsonCurrent.m
%> @brief Example to calculate DC Josephson current through a one-dimensional chain.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @reference
%> E. Perfetto, G. Stefanucci, and M. Cini, Equilibrium and time-dependent Josephson current in one-dimensional superconducting junctions, Phys. Rev. B 80, 205408 (2009)
%> @Available
%> EQuUs v4.5 or later
%> @expectedresult
%> @image html JosephsonCurrent.jpg
%> @image latex JosephsonCurrent.jpg
%> @}
%
%> @brief Example to calculate DC Josephson current through a one-dimensional chain.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function JosephsonCurrent( filenum )

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

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

    % filename containing the input XML    
    inputXML = 'Input_PRB_80_205408.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   = [];    

    % The Fermi energy in the calculations
    EF = 1e-6; % in eV
    % number of energy points on the integration contour
    Edb = 1111;
    
    % The width of the 1D chain
    width = 1;
    % The length of the 1D chain
    height = 8; 
    
    % number of phase difference points on the CPR curve
    phi_db = 40;
    % spacing between the phase difference points
    deltaPhi = pi/(phi_db+1);     
    % 1D array containing the superconducting phase difference values
    DeltaPhi_vec = 0:deltaPhi:pi;      
    %DeltaPhi_vec = [pi/3];  
        
    % *** define the variables storing the output values ***
    
    % results calculated by the scattering technique
     currentvec_continuum_scattering = [];
    current_surf = [];    
    % results calculated by the complex integration method (with Ribbon class)
    currentvec_continuum = [];    
    currentvec_ABS = [];
    currentvec_NBS = [];
    currentvec_1range = [];    
    % results calculated by the complex integration method (with NTerminal class)
    currentvec_TwoTerminal = [];
    
    % creating output directories
    setOutputDir();       
    
    % Calculates the current using the predefined square lattice framework (Ribbon class)
    currentVSDeltaPhi();
    
    % Calculates the current using custom Hamiltonians (TwoTerminal class and CustomHamiltonians class)
    currentVSDeltaPhi_TwoTerminal()
    
    
%% currentVSDeltaPhi
%> @brief Calculates the current using the predefined square lattice framework (Ribbon class)
    function currentVSDeltaPhi()
        
    % setting the phase difference 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_vec(1));
        
    % creating the Ribbon structure describing the junction    
    cRibbon = Ribbon( 'EF', EF, 'width', width, 'height', height, 'silent', true, 'Opt', Opt, 'param', param, 'filenameOut', fullfile(outputdir, [outfilename,'_twoTerminal.xml']));
           
    SNSJosephson_handles = SNSJosephson( Opt, 'junction', cRibbon, 'gfininvfromHamiltonian', false);
    % continuum contribution
    
    % calculate the contribution of the scattering states via the scattering method
    if isempty(currentvec_continuum_scattering)
        [currentvec_continuum_scattering, current_surf] = SNSJosephson_handles.CurrentCalc_continuum( DeltaPhi_vec, 'Edb', Edb );
        save( fullfile(outputdir,[outfilename,'.mat']));
        PlotFunction();
    end
    
    if isempty(currentvec_continuum)
        currentvec_continuum = SNSJosephson_handles.CurrentCalc_discrete( DeltaPhi_vec, 'range', 'CONT', 'Edb', Edb );
        save( fullfile(outputdir,[outfilename,'.mat']));
        PlotFunction();
    end   
    
    % contribution from Andreev bound states
    if isempty(currentvec_ABS)
        currentvec_ABS  = SNSJosephson_handles.CurrentCalc_discrete( DeltaPhi_vec, 'range', 'ABS', 'Edb', Edb );
        save( fullfile(outputdir,[outfilename,'.mat']));
        PlotFunction();
    end

    % contribution from normal bound states
    if isempty(currentvec_NBS)
        currentvec_NBS  = SNSJosephson_handles.CurrentCalc_discrete( DeltaPhi_vec, 'range', 'NBS', 'Edb', Edb );
        save( fullfile(outputdir,[outfilename,'.mat']));
        PlotFunction();
    end
    
    % all the contributions calculated at once
    if isempty(currentvec_1range)
        currentvec_1range = SNSJosephson_handles.CurrentCalc_discrete( DeltaPhi_vec, 'range', 'ALL', 'Edb', Edb );
        save( fullfile(outputdir,[outfilename,'.mat']));
        PlotFunction();
    end
    
    
 
    
    end


%% currentVSDeltaPhi_TwoTerminal
%> @brief Calculates the current using custom Hamiltonians (see the documentation of classes #NTerminal and #Custom_Hamiltonians)
    function currentVSDeltaPhi_TwoTerminal()
        
        % setting the phase difference 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_vec(1));        
        
        % set the parameters to use external source for Hamiltonians
        Opt.Lattice_Type = [];
        Opt.custom_Hamiltonians = 'Custom';
        
        % create an instance of class NTerminal
        cNTerminal = NTerminal('Opt', Opt, 'param', param, 'filenameOut', fullfile(outputdir, [outfilename,'_twoTerminal.xml']), 'WorkingDir',  directory, 'EF', EF, 'CustomHamiltoniansHandle', @Hamiltonians );                   
        
        % create an instance of class SNSJosephson to calculate the DC Josephson current
        cSNSJosephson = SNSJosephson( Opt, 'gfininvfromHamiltonian', true,  'junction', cNTerminal);       

        % perform the calculations
        currentvec_TwoTerminal = cSNSJosephson.CurrentCalc_discrete(DeltaPhi_vec, 'Edb', Edb, 'range', 'ALL');
        
        % plot the results
        PlotFunction();
        
    end


%% Hamiltonians
%> @brief Function to create the custom Hamiltonians for the 1D chain
%> @param varargin Cell array of optional parameters identical to #Custom_Hamiltonians.LoadHamiltonians.
%> @return [1] Cell array of Hamiltonians of one slab in the leads
%> @return [2] Cell array of overlap integrals of one slab in the leads
%> @return [3] Cell array of couplings between the slabs of the leads
%> @return [4] Cell array of overlap integrals between the slabs of the leads
%> @return [5] Cell array of transverse coupling between the unit cells of the lead
%> @return [6] Cell array of #coordinates of the leads
%> @return [7] Hamiltonian of the scattering region
%> @return [8] Overlap integrals of the scattering region
%> @return [9] Transverse coupling for the scattering region
%> @return [10] Overlap integrals for the transverse coupling in the scattering region
%> @return [11] #coordinates of the scattering region
%> @return [12] Cell array of couplings between the scattering region and the leads
%> @return [13] Cell array of overlap integrals between the scattering region and the leads
    function [H0, S0, H1, S1, H1_transverse, coordinates, Hscatter, Sscatter, Hscatter_transverse, Sscatter_transverse, coordinates_scatter, Hcoupling, Scoupling] = Hamiltonians( varargin )
        p = inputParser;
        p.addParameter('q', []);
        p.parse(varargin{:});
        
        q = p.Results.q;
        
        S0        = cell(2,1);
        S1        = cell(2,1);
        Sscatter  = [];
        Scoupling = cell(2,1);
        
        coordinates = cell(2,1);
        coordinates{1} = structures('coordinates');
        coordinates{2} = structures('coordinates');
        coordinates_scatter = structures('coordinates');
        
        [Opt, param] = parseInput('Input_PRB_80_205408.xml');     
        
        H0 = cell(2,1);
        H0{1} = param.Leads{1}.epsilon;
        H0{2} = param.Leads{2}.epsilon;
        
        H1 = cell(2,1);
        H1{1} = param.Leads{1}.vargamma;
        H1{2} = param.Leads{2}.vargamma; 
        
        H1_transverse = cell(2,1);
        H1_transverse{1} = [];
        H1_transverse{2} = [];   
        
        Hscatter_transverse = [];
        Sscatter_transverse = [];
        

        dim = height;
        Hscatter = sparse(1:dim, 1:dim, param.scatter.epsilon, dim, dim) + ...
            sparse(1:dim-1, 2:dim, param.scatter.vargamma, dim, dim) + ...
            sparse(2:dim, 1:dim-1, param.scatter.vargamma, dim, dim);       
        
        Hcoupling = cell(2,1);
        Hcoupling{1} = sparse(1, 1, param.Leads{1}.vargamma, 1, dim);
        Hcoupling{2} = sparse(1, dim, param.Leads{2}.vargamma, 1, dim);
        
                
    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;
          
        
        % creating axes of the plot         
        axes1 = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,... 
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        % Create xlabel
        xlabel('\Delta\phi/\pi','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);
        
        % Create ylabel
        ylabel('I','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);
        
        % define the variable for the legends
        legend_labels = [];
        
        % plot the data
        if ~isempty( currentvec_continuum )
            indexes = logical( currentvec_continuum );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_continuum(indexes), 'b-', 'Parent', axes1)
               legend_labels{end+1} = 'CONT';
        end

        if ~isempty( currentvec_continuum_scattering )
            indexes = logical( currentvec_continuum_scattering );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_continuum_scattering(indexes), 'b--', 'Parent', axes1)
               legend_labels{end+1} = 'CONT S-matrix';
        end
        
        if ~isempty( currentvec_ABS )
            indexes = ~isnan( currentvec_ABS );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_ABS(indexes), 'Linestyle', '-', 'Color', [0 0.83 0], 'Parent', axes1) 
               legend_labels{end+1} = 'Andreev bound states';
        end
        
        if ~isempty( currentvec_NBS )
            indexes = ~isnan( currentvec_NBS );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_NBS(indexes), 'r-', 'Parent', axes1)
               legend_labels{end+1} = 'Normal bound states';
        end
        
        if (~isempty( currentvec_NBS )) && (~isempty( currentvec_continuum )) && (~isempty( currentvec_NBS ))
            total = currentvec_continuum+currentvec_ABS+currentvec_NBS;
            indexes = ~isnan( total );
            plot(DeltaPhi_vec(indexes)/pi, total, 'k-', 'Parent', axes1)
               legend_labels{end+1} = 'Total current';
        end          
        
        
        if ~isempty( currentvec_1range )
            indexes = ~isnan( currentvec_1range );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_1range(indexes), 'c-', 'Parent', axes1)
               legend_labels{end+1} = 'ALL';
        end
        
        
        if ~isempty( currentvec_TwoTerminal )
            indexes = ~isnan( currentvec_TwoTerminal );
            plot(DeltaPhi_vec(indexes)/pi, currentvec_TwoTerminal(indexes), 'Color', [0.8 0 0.6], 'Parent', axes1)
               legend_labels{end+1} = 'TwoTerminal';
        end 
      
           
        % setting the legends
        fig_legend = legend(axes1, legend_labels);
        set(fig_legend, 'FontSize', fontsize, 'FontName','Times New Roman', 'Box', 'off', 'Location', 'NorthEast')

        %set the position of the axis
        figure_pos = get( figure1, 'Position' );
        OuterPosition = get(axes1, 'OuterPosition');
        OuterPosition(1) = 0;
        OuterPosition(2) = 0;
        OuterPosition(3) = figure_pos(3);
        OuterPosition(4) = figure_pos(4);
        set(axes1, 'OuterPosition', OuterPosition);           
                
        
        
        % 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)]; 
        end
               

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



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

        
        
    
end