SNSJosephson.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - SNSJosephson
%    Copyright (C) 2009-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 utilities Utilities
%> @{
%> @file SNSJosephson.m
%> @brief A class to calculate the DC Josephson current.
%> @} 
%> @brief A class to calculate the DC Josephson current.
%%
classdef SNSJosephson < UtilsBase

    properties
          %> The magnetic field in Tesla (OBSOLETE)
        B = 10e-10; %in Tesla      
    end



methods (Access=public)
    
%% Contructor of the class
%> @brief Constructor of the class.
%> @param Opt An instance of the structure #Opt.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
    function obj = SNSJosephson( Opt, varargin )     
        
        obj = obj@UtilsBase( Opt, varargin{:} ); 
        
        obj.T = 0;
        obj.mu = 0;
        obj.scatterPotential = [];
        obj.gfininvfromHamiltonian = false;
        obj.junction = [];
        obj.useSelfEnergy = true;
                
        if strcmpi(class(obj), 'SNSJosephson') 
            obj.InputParsing( varargin{:} );
        end
        
        
    end



%% CurrentCalc_discrete
%> @brief Calculates the Josephson current using the method in PRB 93, 224510 (2016)
%> @param DeltaPhi_vec Array of phase diffrencies between the superconducting contacts.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'range' String. Set 'ABS' (default) for the Andreev bound states, 'NBS' for the normal bound states, 'CONT' for the scattering states, or 'ALL' for calculating all ranges at once. 
%> @param 'Edb' The number of the energy points over the contour integral. (default value is 511)
%> @param 'DeltaPhi' A parameter to control the incident angle of the integration contour near the real axis. (Default value is \f$\Delta\Phi=0.5\pi\f$).
%> @param 'Evec' An array containing the complex energy points in case of custom integration path.
%> @param 'Emin' A minimum of the energy array on the real axis.
%> @param 'Emax' A maximum of the energy array on the real axis.
%> @param 'ProximityOn' Logical value. Set true to recalculate the Green operator of the scattering center for every phase difference. Useful for self-consistent modelling of the proximity effect.
%> @return An array of the calculated current corresponding to the phase differencies DeltaPhi_vec
    function currentvec = CurrentCalc_discrete( obj, DeltaPhi_vec, varargin )
        
        p = inputParser;
        p.addParameter('range', 'ABS' ); % NBS: normal bound states, ABS: andreev bound states, CONT: scattering range, ALL: all contributions at once
          p.addParameter('Edb', 511 ); % number of energy points on the contour 
        p.addParameter('DeltaPhi', 0.5*pi );
        p.addParameter('Evec', [] );
        p.addParameter('Emin', [] );
        p.addParameter('Emax', [] );
        p.addParameter('ProximityOn', false ); %New if proxymity effect is taken into account, the current operator must be recalculated for each phase difference
        
        p.parse(varargin{:});       
        range   = p.Results.range;
        Edb      = p.Results.Edb;
        Evec     = p.Results.Evec;
        Emin     = p.Results.Emin;
        Emax     = p.Results.Emax;
        DeltaPhi = p.Results.DeltaPhi;
        ProximityOn = p.Results.ProximityOn;
        
        DeltaPhi_vec = [0, DeltaPhi_vec];
        DeltaPhi_vec = diff(DeltaPhi_vec);
        
        % reset the phase difference to zero
        obj.ResetPhase()
        

        if ProximityOn
            obj.gfininvfromHamiltonian = true;
        end
        
        
        phivec = (0.5 + 0.5*sin( -0.5*pi:pi/(Edb+1):0.5*pi )) * DeltaPhi + (pi-DeltaPhi)/2;
        
        if strcmp( range, 'NBS' ) && isempty(Evec) && (isempty(Emax) && isempty(Emin))
            obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Calculating current for the normal bound states.'], 1)
            obj.getBandWidth();
            Emax = -obj.BandWidth.Lead.Emax;
            Emin = -obj.BandWidth.Scatter.Emax;
%            phivec = (0.5 + 0.5*sin( -0.5*pi:pi/(Edb+1):0.5*pi )) * pi;
            %fact = -fact; %the contour integral is in the opposite direction as for the ABS
        elseif strcmp( range, 'ABS' ) && isempty(Evec) && (isempty(Emax) && isempty(Emin))
            obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Calculating current for the Andreev bound states.'], 1)
            obj.getBandWidth();
            param = obj.junction.FL_handles.Read( 'param' );         
            pair_potential = min([abs(param.Leads{1}.pair_potential), abs(param.Leads{2}.pair_potential)]);
            Emin = -pair_potential;
            Emax = -0;
            %phivec = (0.5 + 0.5*sin( -0.5*pi:pi/(Edb+1):0.5*pi )) * pi;
        elseif strcmp( range, 'CONT' ) && isempty(Evec) && (isempty(Emax) && isempty(Emin))
            obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Calculating current for the scattering states.'], 1)
            obj.getBandWidth();
            param = obj.junction.FL_handles.Read( 'param' );         
            pair_potential = min([abs(param.Leads{1}.pair_potential), abs(param.Leads{2}.pair_potential)]);
            Emin = -obj.BandWidth.Lead.Emax;
            Emax = -pair_potential;
            %phivec = (0.5 + 0.5*sin( -0.5*pi:pi/(Edb+1):0.5*pi )) * pi;    
        elseif strcmp( range, 'ALL' ) && isempty(Evec) && (isempty(Emax) && isempty(Emin))
            obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Calculating current for all the states.'], 1)
            obj.getBandWidth();        
            Emin = min([-obj.BandWidth.Lead.Emax,-obj.BandWidth.Scatter.Emax]);
            Emax = -0;
            %phivec = (0.5 + 0.5*sin( -0.5*pi:pi/(Edb+1):0.5*pi )) * pi;        
        elseif ~isempty(Evec)
            obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Calculating current for custom energy contour.'], 1)  
            Emin = min(real(Evec));
            Emax = max(real(Evec)); 
            Edb = numel(Evec);
        elseif (~isempty(Emax) && ~isempty(Emin))
        else
            err = MException(['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete', 'Unrecognized energy range']);
            save('Error_SNSJosephson_CurrentCalc_discrete.mat');
            throw(err); 
        end
        
        if Emin >= Emax
            currentvec = zeros(size(DeltaPhi_vec));
            return
        end    
        
        % E -> -E symmetry for both zero and nonzero temperature.
        fact = 2;
               
        % creating the energy array
        radius = abs( Emax-Emin) /2;  
        R = radius/sin(DeltaPhi/2);
        if isempty(Evec)
            Evec = R*(cos(phivec) + 1i*sin(phivec)) - radius + Emax - 1i*R*sin((pi-DeltaPhi)/2);
        end
        
        obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Calculating the current between Emin = ', num2str(Emin),' and Emax = ',num2str(Emax),],1);
          deltaEvec = [0, diff(Evec)];
        currentsurf = complex(zeros( length(Evec), length(obj.T), length(DeltaPhi_vec) ));
        
        % open the parallel pool
        poolmanager = Parallel( obj.Opt );
        poolmanager.openPool(); 

        % setting function handles for the parallel calculations
        hdisplay = @obj.display;
        hSetEnergy = @obj.SetEnergy;
        hcreateCurrentOperator = @obj.createCurrentOperator;
        hincreasePhaseDifference = @obj.increasePhaseDifference;
        hcomplexCurrent = @complexCurrent;
        
        junction_loc = obj.junction;
        T = obj.T;
        
        % starting parallel calculations
        parfor (iidx = 1:length(Evec), obj.Opt.workers)
        %for iidx = 1:length(Evec)
            
if isempty(junction_loc.getTransverseMomentum())            
            tic                          
            feval( hdisplay, ['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Currently evaluating ', num2str(iidx), ' from ', num2str(length(Evec))], 1);
end
            
            % in the first iteration the self energies of all terminals should be recalculated
               recalculateSurface = [1 2];
            
            % set the energy value in the current iteration
            feval(hSetEnergy, Evec(iidx), 'junction', junction_loc );
            
            try
                % determine the current operator and calculate the Green operator of the scattering center
                current_op = feval(hcreateCurrentOperator, 'junction', junction_loc);
            catch errCause
                currentsurf( iidx, :, : ) = NaN(length(obj.T), length(DeltaPhi_vec));
                feval(hdisplay, ['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Error in SNSJosephson:CurrentCalc_discrete at energy point E = ', num2str(Evec(iidx)), ' caused by:'],1);
                feval(hdisplay, errCause.identifier, 1 );
                    feval(hdisplay, errCause.stack(1), 1 );
                feval(hdisplay, ['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Giving NaN for the current at the given energy point'], 1); 
               continue
            end 
            
            % preallocating array for the calculated, phase resolved current values
            currentvec2int = zeros( length( T ), length( DeltaPhi_vec ) );
            for jjdx = 1:length( DeltaPhi_vec )                
                try 
                    % increasing the phase different between the superconducting leads
                    feval(hincreasePhaseDifference, DeltaPhi_vec(jjdx), junction_loc );
                    if ProximityOn && jjdx > 1
                         % if proximity effect is taken into account, the current operator must be always recalculated
                         current_op = feval(hcreateCurrentOperator, 'junction', junction_loc);
                    end
                    
                    % evaluate the current at a complex energy
                    currentvec2int( :, jjdx ) = feval(hcomplexCurrent, junction_loc, current_op, recalculateSurface);
                    
                    % the phase different is increased by unitary transformation, so recalculating the leeds is not necessary
                    recalculateSurface = []; % see increasePhaseDifference
                catch errCause
                    feval(hdisplay, ['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Error in SNSJosephson:CurrentCalc_discrete at energy point E = ', num2str(Evec(iidx)), 'and at phase difference ', num2str(DeltaPhi_vec(jjdx)), ' caused by:'], 1);
                    feval(hdisplay, errCause.identifier, 1 );
                    for jjjdx = 1:length(errCause.stack)
                        feval(hdisplay, errCause.stack(jjjdx), 1 );
                    end
                    feval(hdisplay, ['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Giving NaN for the calculated current at the given energy point and phase difference.'], 1); 
                    currentvec2int( jjdx ) = NaN;
                end
                
            end 
            currentsurf( iidx, :, : ) = currentvec2int;
            
if isempty(junction_loc.getTransverseMomentum())                  
            toc
end
        end

        % closing the parallel pool
        poolmanager.closePool();
          
        
          % calculating the contour integrals
        currentvec = cell( size(obj.T) );
          for kkdx = 1:length( obj.T )
               currentvec{kkdx} = zeros(size(DeltaPhi_vec));          
             for jjdx = 1:length( DeltaPhi_vec )
                    currentvec2integrate = currentsurf(:,kkdx,jjdx);
                    % omitting NaNs by substituting interpolates
               indexes = isnan(currentvec2integrate); 
                    indexesall = 1:length(currentvec2integrate);
                       
               currentvec2integrate(indexes) = interp1(indexesall(~indexes), currentvec2integrate(~indexes), indexesall(indexes), 'spline' );
                    current = fact*imag( obj.SimphsonInt( currentvec2integrate.*transpose(deltaEvec) ))/(pi);
                    currentvec{kkdx}( jjdx ) = currentvec{kkdx}( jjdx ) + current;
               end
        end    
          if length(obj.T) == 1
               currentvec = currentvec{1};
          end
        
        obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: Process finished.'], 1)
        

        % nested functions


          %------------------------------------------------
          % calculates current for complex energy
          function current = complexCurrent( junction_loc, current_op, recalculateSurface )
            
        
            if isprop(junction_loc, 'Ribbons')
                junction_loc = junction_loc.Ribbons{1}; %for the CombineRibbons interface
            end                    
                              
            [~, gfininv  ] = junction_loc.GetFiniteGreensFunction();
            G = junction_loc.CustomDysonFunc('recalculateSurface', recalculateSurface, ...
                'decimate', true, 'constant_channels', false, ...
                'SelfEnergy', obj.useSelfEnergy, 'keep_sites', 'scatter');
             
               G_size = size(G,1);                               
               gfininv_size = size(gfininv,1);                   
               if G_size ~= gfininv_size
                    warning( ['EQuUs:Utils:', class(obj), ':CurrentCalc_discrete: wrong size of Greens function'] )
                    current = 0;
                    return
               end

               fE = obj.Fermi(junction_loc.getEnergy());
             
            scatter_sites = junction_loc.CreateH.Read('kulso_szabfokok');            
            if ~isempty(scatter_sites) 
                % if the Green operator is calculated from the whole Hamiltonian, the coordinates contain all sites.
                coordinates = junction_loc.CreateH.Read('coordinates');
                BdG_indexes = coordinates.BdG_u(scatter_sites);
            else
                % obtaining only the surface points of the junction structure
                coordinates = junction_loc.getCoordinates();
                BdG_indexes = coordinates.BdG_u();
            end
            
               current = zeros(length( fE ), 1 );

            for idx = 1:length( fE )
                    G_tmp = G;
                    G_tmp( BdG_indexes, BdG_indexes) = G_tmp( BdG_indexes, BdG_indexes)*fE(idx);
                    G_tmp( ~BdG_indexes, ~BdG_indexes) = G_tmp( ~BdG_indexes, ~BdG_indexes)*(1-fE(idx));

                    current(idx) = trace( current_op*G_tmp );

            end
         
        end
        
        % end nested functions
  
    end


%% CurrentCalc_continuum
%> @brief Calculates the Josephson current of the continous scattering states.
%> @param DeltaPhi_vec Array of phase diffrencies between the superconducting contacts.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'Edb' The number of the energy points over the contour integral. (default value is 511)
%> @return [1] An array of the calculated current corresponding to the phase differencies DeltaPhi_vec
%> @return [2] Energy resolved Josephson current.
    function [currentvec, current_surf] = CurrentCalc_continuum( obj, DeltaPhi_vec, varargin  )
        
        obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_continuum: Calculating current in the continuum range from the scattering matrix.'], 1)
        p = inputParser;
          p.addParameter('Edb', 511 ); % number of energy points on the contour 
        p.parse(varargin{:});       
          Edb     = p.Results.Edb;

        DeltaPhi_vec = [0, DeltaPhi_vec];
        DeltaPhi_vec = diff(DeltaPhi_vec);        
     currentvec = zeros(size(DeltaPhi_vec));
        
        % reset the phase difference to zero
        obj.ResetPhase()
        
     % determine the bandwidth
     obj.getBandWidth();
        Emax = -obj.BandWidth.Lead.Emin;
        Emin = -obj.BandWidth.Lead.Emax;
        Evec = (0.5+0.5*sin((0:Edb-1)/(Edb-1)*pi - pi/2))*(Emax-Emin) + Emin;
        deltaEvec = [0, diff(Evec) ];
        
        current_surf = zeros(length(DeltaPhi_vec), length(Evec));
        % E -> -E symmetry for both zero and nonzero temperature
        fact = 2;
        
        hCurrentCalc_continuum_forGivenEnergy = @CurrentCalc_continuum_forGivenEnergy;
        ribbon_loc = obj.junction;%;.CreateClone();    
        for kdx = 1:length( Evec )
            obj.display( ['EQuUs:Utils:', class(obj), ':CurrentCalc_continuum: Evaluating ', num2str(kdx), ' from ', num2str(length(Evec))], 1);
tic            
            current_surf(:,kdx) = feval(hCurrentCalc_continuum_forGivenEnergy, DeltaPhi_vec, Evec(kdx), ribbon_loc ); 
toc
        end
        
        for jdx = 1:length( DeltaPhi_vec )
            currentvec2integrate =  current_surf(jdx,:);
            % omitting NaNs by substituting interpolates
          indexes = isnan(currentvec2integrate); 
               indexesall = 1:length(currentvec2integrate);
                       
            currentvec2integrate(indexes) = interp1(indexesall(~indexes), currentvec2integrate(~indexes), indexesall(indexes), 'spline' );
            current_surf(jdx,:) = currentvec2integrate;
               currentvec(jdx) = currentvec(jdx) + fact*obj.SimphsonInt( currentvec2integrate.*(deltaEvec) )/(2*pi);            
        end
        
        %reset the system  for a new calculation
        obj.InputParsing();
        
        obj.display(['EQuUs:Utils:', class(obj), ':CurrentCalc_continuum: Process finished.'], 1)
        
        %-----------------------------------------------------------
        % nested functions
        function current = CurrentCalc_continuum_forGivenEnergy( DeltaPhi_vec, E, junction_loc )
        
            obj.SetEnergy( E, 'junction', junction_loc  );  
        
            param = junction_loc.FL_handles.Read( 'param' );         
            pair_potential = min([abs(param.Leads{1}.pair_potential), abs(param.Leads{2}.pair_potential)]);
        
            current = NaN( size( DeltaPhi_vec ) );
            if abs(obj.junction.getEnergy()) <= abs(pair_potential)
                return
            end
        
            current_op = obj.createCurrentOperator('junction', junction_loc); 
            if isnan(current_op)
                return
            end      
              
            recalculateSurface = [1 2];
            for idx = 1: length( DeltaPhi_vec )
                try   
                
                    % setting phi2 to phi2+delta_phi 
                    obj.increasePhaseDifference( DeltaPhi_vec(idx), junction_loc );
                    %calculating the current
                    current_tmp = CurrentCalc_continuum_for_GivenDeltaPhi( junction_loc, current_op, recalculateSurface);
                catch errCause
                    obj.display( ['EQuUs:Utils:', class(obj), ':CurrentCalc_continuum: Error in SNSJosephson:CurrentCalc_continuum_forGivenEnergy caused by:'])
                    obj.display( errCause.identifier, 1 );
                    obj.display( errCause.stack(1), 1 );
                    obj.display( ['EQuUs:Utils:', class(obj), ':CurrentCalc_continuum: Giving NaN for the calculated current'])
                    current_tmp = NaN;
                end
                recalculateSurface = [];
                current(idx) = current_tmp;
            end       
        
            
    
            %-------------------------------------------------
            function current = CurrentCalc_continuum_for_GivenDeltaPhi( ribbon_loc, current_op, recalculateSurface )
  
                %Dysonfunc = @()CustomDysonFunc('recalculateSurface', recalculateSurface);
                Dysonfunc = @()ribbon_loc.CustomDysonFunc('recalculateSurface', recalculateSurface, ...
                    'decimate', false, 'constant_channels', false, 'SelfEnergy', obj.useSelfEnergy, 'keep_sites', 'scatter' );
            
                G = ribbon_loc.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
                        
                %[S,ny] = obj.junction.FL_handles.SmatrixCalc();
                M = ribbon_loc.FL_handles.Read('M');
            
                % The calculated leads for calculate the scattering states
                Leads = ribbon_loc.FL_handles.Read('Leads');
                vcsop_p = cell(length(Leads),1);
                expk_p = cell(length(Leads),1);
                expk_m = cell(length(Leads),1);
                modusmtx_p = cell(length(Leads),1);
                Normamtx = cell(length(Leads),1);
                modusmtx_m = cell(length(Leads),1);
                d_modusmtx_m = cell(length(Leads),1);
                for iidx = 1:length(Leads)
                    vcsop_p{iidx}    = Leads{iidx}.Read('vcsop_p');
                    expk_p{iidx}        = Leads{iidx}.Read('expk_p');
                    expk_m{iidx}        = Leads{iidx}.Read('expk_m');
                    modusmtx_p{iidx} = Leads{iidx}.Read('modusmtx_p');
                    Normamtx{iidx}   = Leads{iidx}.Read('Normamtx');
                    modusmtx_m{iidx} = Leads{iidx}.Read('modusmtx_m');
                    d_modusmtx_m{iidx} = Leads{iidx}.Read('d_modusmtx_m');
                end
            
                fE = obj.Fermi(ribbon_loc.getEnergy());
          
                %propagating indexes
                tolerance = 1e-6;
                prop_indexes{1} = find(abs(abs(expk_p{1}) - 1) < tolerance);
                prop_indexes{2} = find(abs(abs(expk_p{2}) - 1) < tolerance);

                current = 0;
                for cell_index = 1:length(Leads)
                    for iidx = 1:length(prop_indexes{cell_index})
                        current_tmp = CurrentFromGivenMode( cell_index, prop_indexes{cell_index}(iidx) );
                        current = current + current_tmp;
                    end
                end    
       
                if imag(current) > 1e-6
                    warning(['EQuUs:Utils:', class(obj), ':CurrentCalc_continuum: Current is complex: ', num2str(current)])
                end
                current = real(current);
                
        
                %--------------------------------------------------
               % calculates current from given incomming mode
                function current_trans = CurrentFromGivenMode( cell_index, idx )
                    tolerance = 1e-6;
                    modus_in = zeros(size(G,1), 1);
                    indexes = sum(M(1:cell_index-1)) + (1:M(cell_index));
                    modus_in(indexes) = Normamtx{cell_index}*modusmtx_p{cell_index}(:,idx)/sqrt(vcsop_p{cell_index}(idx));
            
                    modus_out = G*modus_in;
                    sc_wavefnc = modus_out(sum(M)+1:end);   
                
                    coordinates = ribbon_loc.getCoordinates();
                    scatter_sites = ribbon_loc.CreateH.Read('kulso_szabfokok');
                    if ~isempty(scatter_sites) % if the Green operator is calculated from the full Hamiltonian, the coordinates contain all sites.
                        BdG_indexes = coordinates.BdG_u(scatter_sites);
                    else
                        BdG_indexes = coordinates.BdG_u;
                    end
            
                    current_op_particle = zeros(size(current_op));
                    current_op_particle( BdG_indexes, BdG_indexes) = current_op( BdG_indexes, BdG_indexes)*fE;
                    current_op_particle( ~BdG_indexes, ~BdG_indexes) = current_op( ~BdG_indexes, ~BdG_indexes)*(1-fE);  
            
                
%G_tmp = G(sum(M)+1:end, sum(M)+1:end);
%                   G_tmp( BdG_indexes, BdG_indexes) = G_tmp( BdG_indexes, BdG_indexes)*fE;
%                   G_tmp( ~BdG_indexes, ~BdG_indexes) = G_tmp( ~BdG_indexes, ~BdG_indexes)*(1-fE);
                              
%                   current_trans_old = imag(trace( current_op*G_tmp )); 
                
%G_tmp = G(sum(M)+1:end, sum(M)+1:end);


                %current_trans = imag(trace( current_op_particle*G_tmp ));
                %return
            
%               particle current
                    current_trans = sc_wavefnc'*current_op_particle*sc_wavefnc;
                    return

                    %transmitted quasiparticle current      
                    quasi_current_trans = sc_wavefnc'*current_op*sc_wavefnc;
                
                    %H1sc = ribbon_loc.Scatter_UC.Read('H1');
                    %c1 = sc_wavefnc(1:M(cell_index));
                    %c2 = sc_wavefnc(M(cell_index)+1:2*M(cell_index));
                    %quasi_current_trans = 1/(1i)*( c2'*H1sc'*c1 - c1'*H1sc*c2)
            
                    % ******** unitarity test **********
                    % coupling Hamiltonian to calculate the current in the incomming lead
                    H1 = Leads{cell_index}.Read('H1');     

                    modus_in = modusmtx_p{cell_index}(:,idx)/sqrt(vcsop_p{cell_index}(idx));
                    quasi_current_in = 1i*modus_in'*(H1*expk_p{cell_index}(idx)-H1'/expk_p{cell_index}(idx))*modus_in;
            
                    modus_ref = modus_out(indexes) - modus_in;
                    r = d_modusmtx_m{cell_index}*modus_ref;
                    quasi_current_ref = 0;
                    for iiidx = 1:M(cell_index)
                        if abs(abs(expk_m{cell_index}(iiidx))-1) <= 1e-6
                            modus_ref_tmp = r(iiidx)*modusmtx_m{cell_index}(:,iiidx);
                            quasi_current_ref = quasi_current_ref + 1i*modus_ref_tmp'*(H1*expk_m{cell_index}(iiidx)-H1'/expk_m{cell_index}(iiidx))*modus_ref_tmp;
                        end
                    end
                    quasi_current_ref = -quasi_current_ref;
            
                    delta_quasi_current = real(quasi_current_in - quasi_current_trans - quasi_current_ref);
                    delta_quasi_current
                    if abs(delta_quasi_current) > tolerance
                        display(['EQuUs:Utils:', class(obj), ':CurrentCalc_continuum: unitarity test failed with \Delta I = ',num2str(delta_quasi_current)]);
                    end 
                    % ******** unitarity test **********
            
                end
        
            end
    
        end
        % end nested functions
    end

%% createCurrentOperator
%> @brief Calculates the charge current operator from the inverse of the Greens function of the scattering region.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'junction' An instance of class #NTerminal (or its subclass) describing the junction.
%> @return [1] The matrix operator of the current operator.
     function current_op = createCurrentOperator( obj, varargin )
        
        p = inputParser;
        p.addParameter('junction', obj.junction);
        p.parse(varargin{:});
        
        junction_loc = p.Results.junction;
          
          obj.create_scatter_GreensFunction('junction', junction_loc, 'gauge_trans', true);
        [~, gfininv  ] = junction_loc.GetFiniteGreensFunction();
        if isnan(gfininv)
               current_op = NaN;
               return
        end 
        
        if strcmpi( class(junction_loc), 'Ribbon' )
            M = size(gfininv)/4;
     
            current_op = zeros( size(gfininv) );
            % gfininv = E - Heff(E) !!!!!
            current_op( 1:M, 2*M+1:3*M ) = -1/1i*(-1)*gfininv( 1:M, 2*M+1:3*M );
            current_op( M+1:2*M, 3*M+1:end ) = -1/1i*(-1)*gfininv( M+1:2*M, 3*M+1:end );
            current_op( 2*M+1:3*M, 1:M ) = 1/1i*(-1)*gfininv( 2*M+1:3*M, 1:M );
            current_op( 3*M+1:end, M+1:2*M ) = 1/1i*(-1)*gfininv(3*M+1:end, M+1:2*M );        
        else
        % gfininv = E - Heff(E) !!!!!
            current_op = 1/1i*(triu(gfininv) - tril(gfininv)); 
        end

    end


%% getBandWidth
%> @brief Determines the band width of the leads and of the scattering region.
%> @return An instance of structure #BandWidth containing the bandwidths.
    function ret = getBandWidth( obj )
     obj.display('Getting Bandwidths', 1);
        
        if ~isempty(obj.BandWidth)
            return
        end
        
        poolmanager = Parallel( obj.Opt );
        poolmanager.openPool();
        
        Surface_1 = obj.junction.FL_handles.SurfaceGreenFunctionCalculator(1, 'Just_Create_Hamiltonians', true, 'q', obj.junction.getTransverseMomentum(), 'leadmodel', obj.junction.leadmodel, ...
            'CustomHamiltonian', obj.junction.cCustom_Hamiltonians);
        W = BandWidthCalculator( Surface_1 );
        
        param = obj.junction.FL_handles.Read( 'param' );         
        pair_potential = min([abs(param.Leads{1}.pair_potential), abs(param.Leads{2}.pair_potential)]);
        
        obj.BandWidth.Lead = structures('BandWidthLimits');
        obj.BandWidth.Lead.Emin = abs(pair_potential);
        obj.BandWidth.Lead.Emax = W;
        
        
        if strcmpi( class(obj.junction), 'Ribbon' )
            obj.junction.CreateRibbon();        
            W = BandWidthCalculator( obj.junction.Scatter_UC );
        elseif strcmpi( class(obj.junction), 'NTerminal' )
            obj.junction.CreateScatter();
            W = 1.3*BandWidthCalculatorMolecule( obj.junction.CreateH ); % multiplied by 1.3 just for sure to include normal bound states
        else
            W = 0;
        end
        obj.BandWidth.Scatter = structures('BandWidthLimits');
        obj.BandWidth.Scatter.Emin = 0;
        obj.BandWidth.Scatter.Emax = W;
        
        poolmanager.closePool();

        ret = obj.BandWidth;
         
        % nested functions
         %--------------------------------------------
         function W = BandWidthCalculator( Scatter_UC )

          H0 = Scatter_UC.Read('H0');
            q = obj.junction.getTransverseMomentum();
            if ~isempty( q )
                H1_transverse = Scatter_UC.Read('H1_transverse');
                H0 = H0 + H1_transverse*exp(1i*q) + H1_transverse'*exp(-1i*q); 
            end
          H1 = Scatter_UC.Read('H1'); 
          u = rand(size(H0,1)); 
          u = u/norm(u);     
            
            hgetMaxEigenvalue = @getMaxEigenvalue;
            hsecular_H = @secular_H;
            k = -pi:2*pi/50:pi;
            arnoldi = zeros( length(k),1);
            parfor (idx = 1:length(k), poolmanager.NumWorkers)
                Heff = feval(hsecular_H,H0,H1,k(idx));
                arnoldi(idx) = abs(feval(hgetMaxEigenvalue, Heff ));                
            end
            W = max( arnoldi );
            

          % --------------------------------------
          % Hamiltonian for the secular equation
          function H = secular_H(H0,H1,k)     
               H = H0 + H1*exp(1i*k) + H1'*exp(-1i*k);    
          end
         
            %--------------------------------------------
          % getting the maximal eigenvalue with Arnoldi iteration
          % Numerical Linear Algebra" by Trefethen and Bau. (p250) or
          % http://en.wikipedia.org/wiki/Arnoldi_iteration
          function ret = getMaxEigenvalue( Heff )
               %u = rand(size(H0,1)); 
               %u = u/norm(u);
               maxIter = 1e2;
               tolerance = 1e-5;
              ret = 0;
               lambda = 1;
               for iidx=1:maxIter
                    u_new = Heff*u;
                    norm_new = norm(u_new);
                    lambda_new = norm_new/norm(u);
               
                    if abs((lambda - lambda_new)/lambda) < tolerance
                         ret = lambda;
                         return
                    end                 
                    u = u_new/norm_new;
                    lambda = lambda_new;          
               end
          end
         end
         
         
         %--------------------------------------------
         function W = BandWidthCalculatorMolecule( CreateH )
            Hscatter  = CreateH.Read('Hscatter');
            W = abs(eigs(Hscatter,1));
         end   
         
         % end nested functions
        
    end
    
    
    end % public methods
    
    
    methods (Access=protected)    




%% create_scatter_GreensFunction
%> @brief Calculates the surface Green operator of the scattering region.
%> The calculated Green operator is stored within the objec #junction (or in the object junction given as an optional parameter)
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'gauge_trans' Logical value. Set true (default) to perform gauge transformation on the Green's function and on the Hamiltonians, or false otherwise.
%> @param 'junction' An instance of class #NTerminal (or its subclass) describing the junction.
    function create_scatter_GreensFunction(  obj, varargin )
        
        p = inputParser;
        p.addParameter('gauge_trans', true); % logical: true if want to perform gauge transformation on the Green's function and Hamiltonians
        p.addParameter('junction', obj.junction); %for parallel computations
        p.parse(varargin{:});       
        gauge_trans     = p.Results.gauge_trans;
        ribbon_loc     = p.Results.junction;       
        
        if obj.gfininvfromHamiltonian
            %calculating the Greens function from the Hamiltonian
            ribbon_loc.CalcFiniteGreensFunctionFromHamiltonian('onlyGinv', true, 'gauge_trans', gauge_trans, 'PotInScatter', obj.scatterPotential );
        else
            % Calculating the Greens function by the fast way optimized for translational invariant systems
            ribbon_loc.CalcFiniteGreensFunction('onlyGinv', true, 'gauge_trans', gauge_trans);
        end         
        
        
    end



%% increasePhaseDifference
%> @brief Increase the phase difference between the Leads by delta_phi (apply gauge transformation on the second lead)
%> @param delta_phi The phase difference increment between the leads.
%> @param ribbon_loc An instance of class NTerminal (or its subclass) describing the junction.
    function increasePhaseDifference( obj, delta_phi, junction_loc )
        
            
        if isprop(junction_loc, 'Ribbons')
            junction_loc = junction_loc.Ribbons{1}; %for the CombineRibbons interface
        end
        
        Leads = junction_loc.FL_handles.Read( 'Leads' ); 
        
        if isempty( Leads )
            return
        else
            gaugetrans( Leads{2} ) 
        end    
        
        % gauge transformation in the interface region
        if isempty( junction_loc.Interface_Regions )
            return
        else
            gaugetrans( junction_loc.Interface_Regions{2} ) 
        end   
        
        % nested functions
        
        %--------------------------------------
        % function that performs tha gauge transformation on the pair
        % potential
        function gaugetrans( Lead )
            if abs(delta_phi ) < 1e-16
                return
            end       
                        
            coordinates = Lead.Read('coordinates');
            regular_sites = Lead.Read('kulso_szabfokok');
            
            if isempty(coordinates)
                % Hamiltonians are not created yet  
                return
            end            
            
            if isempty(coordinates.BdG_u)
                % not a BdG calculation
                warning(['EQuUs:Utils:', class(obj), ':gaugetrans:Not a BdG calculation. Gauge transformation of the pair potential is skipped']);
                return
            end
            
            if isempty(regular_sites)
                regular_sites = 1:Lead.Get_Neff();
            end
            
            Umtx_diag = zeros(size(regular_sites));            
        
            Umtx_diag( coordinates.BdG_u(regular_sites) ) = exp(1i*delta_phi/2);
            Umtx_diag( ~coordinates.BdG_u(regular_sites) ) = exp(-1i*delta_phi/2);        
            Umtx = diag(Umtx_diag);
            
            Lead.Unitary_Transform( Umtx );                 
            
        end
        
        % end nested functions
        
    end


%% ResetPhase
%> @brief Resets the phase difference to the initial value.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'junction' An instance of class #NTerminal (or its subclass) describing the junction.
    function ResetPhase( obj, varargin ) 
        p = inputParser;
        p.addParameter('junction', obj.junction);
        p.parse(varargin{:});
                
        junction_loc = p.Results.junction;
        
        % retriving the param structure
        param = junction_loc.getParam();
        
        % setting the phase of the second lead
        junction_loc.setParam( param );
        
    end



    

%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'T' The temperature in Kelvin (scalar or an array) 
%> @param 'mu' The Chemical potential in the same unit as other energy scales in the Hamiltonians.
%> @param 'scatterPotential' A function handle pot=f( #coordinates ) or pot=f( #CreateHamiltonians, Energy) for the potential to be applied in the Hamiltonian (used when FiniteGreensFunctionFromHamiltonian=true).
%> @param 'gfininvfromHamiltonian' Logical value. Set true calculate the surface Greens function of the scattering region from the Hamiltonaian of the scattering region, or false (default) to calculate it by the fast way (see Phys. Rev. B 90, 125428 (2014) for details).
%> @param 'useSelfEnergy' Logical value. Set true to use the self energies of the leads in the Dyson equation, or false (default) to use the surface Green function instead.
%> @param 'junction' An instance of class #NTerminal (or its derived class) representing the junction.
    function InputParsing( obj, varargin )
                
        p = inputParser;
        p.addParameter('T', obj.T);
        p.addParameter('mu', obj.mu, @isscalar);
        p.addParameter('scatterPotential', obj.scatterPotential);
        p.addParameter('gfininvfromHamiltonian', obj.gfininvfromHamiltonian);
        p.addParameter('ribbon', obj.junction);     % NEW OBSOLETE use parameter junction istead
        p.addParameter('junction', obj.junction);     % NEW
        p.addParameter('useSelfEnergy', obj.useSelfEnergy); %NEW
        
        p.parse(obj.varargin{:});
        
        junction             = p.Results.junction;
        if ~isempty(p.Results.ribbon) && isempty(junction)
            junction             = p.Results.ribbon;
        end
            
        
        InputParsing@UtilsBase( obj, 'T', p.Results.T, ...
                                'mu', p.Results.mu, ...
                                'scatterPotential', p.Results.scatterPotential, ...
                                'gfininvfromHamiltonian', p.Results.gfininvfromHamiltonian, ...
                                'useSelfEnergy', p.Results.useSelfEnergy, ...
                                'junction', junction);
        
        
        
        
        obj.junction.FL_handles.LeadCalc('createCore', true, 'q', obj.junction.getTransverseMomentum() );
               

    end 
    
    
end % protected methods
end