test_Keldysh.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities
%    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 unit_tests Unit Tests
%> @{
%> @file test_Keldysh.m
%> @brief Testfile to check functionalities of the classes developed for steady-state Keldysh formalism.
%> @Available
%> EQuUs v4.9 or later
%> @}
%> @brief Testfile to check functionalities of the classes developed for steady-state Keldysh formalism.
function test_Keldysh()

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

    % create input structures
    [Opt, param] = parseInput( 'Input_Ribbon.xml');
    
    % defining the width of the lead 1
    param.Leads{1}.M = 10;
    

    % energy value used in the calculations
    Energy = -0.01;
    
    % The temperature in Kelvins
    T = 10;
    
    %% test of the class Lead_Keldysh
    disp( '***** test of the class Lead_Keldysh *****' );

    % create an instance of class Lead_Keldysh
    cLead_K = Lead_Keldysh(Opt, param, 'Hanyadik_Lead', 1, 'T', T, 'bias', 0.04);

    % create Hamiltonians
    cLead_K.CreateHamiltonians();     
    %cLead_K.CalcSpektrum('ka_min', -1, 'ka_max', 1, 'ka_num', 100, 'toPlot', true)
     
    % Trukkos sajatertek
    cLead_K.TrukkosSajatertekek(Energy);
     
    % Calculate The grou velocity
    cLead_K.Group_Velocity();
       
    % calculate the lesser Self energy
    cLead_K.LesserSelfEnergy();
    
     
    
    %% test of class Ribbon_Keldysh, Transport_Interface_Keldysh
    disp( '***** test of class Ribbon_Keldysh, Transport_Interface_Keldysh *****' );    
    
    % chemical potentials in the leads
    bias_leads = [0.0 0.0];
    
    % chemical potential in the central device
    EF = 0.05;
    
    % create an instance of class Lead_Keldysh
    cRibbon_K = Ribbon_Keldysh('width', 10, 'height', 10, 'Opt', Opt, 'param', param, 'E', Energy, 'bias_leads', bias_leads, 'EF', EF, 'T', T);    
    
    % calculate the inverse of the retarded surface Green operator of the isolated scattering center
    cRibbon_K.CalcFiniteGreensFunctionFromHamiltonian('onlyGinv', true);
    
    % calculate the lesser Green operator
    lesserG = cRibbon_K.LesserGreenFunction();
    
    retardedG = cRibbon_K.FL_handles.Read('G');
    occupation_number = cRibbon_K.Fermi( Energy );
    disp(['Occupation number = ', num2str( occupation_number ), ', Energy: ', num2str(Energy), ', temperature = ', num2str(T)]);
    
    checkpoint = norm( lesserG - occupation_number*(retardedG'-retardedG ) );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
    % set new temperature value
    T = 1000;    
    cRibbon_K.setTemperature( T );
    
    % calculate the lesser Green operator
    lesserG = cRibbon_K.LesserGreenFunction();
    
    retardedG = cRibbon_K.FL_handles.Read('G');
    occupation_number = cRibbon_K.Fermi( Energy );
    disp(['Occupation number = ', num2str( occupation_number ), ', Energy: ', num2str(Energy), ', temperature = ', num2str(T)]);
    
    checkpoint = norm( lesserG - occupation_number*(retardedG'-retardedG ) );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
    
    % set new energy value
    Energy = 0.15;
    
    cRibbon_K.setEnergy( Energy );
    
    % calculate the inverse of the retarded surface Green operator of the isolated scattering center
    cRibbon_K.CalcFiniteGreensFunctionFromHamiltonian('onlyGinv', true);
    
    % calculate the lesser Green operator
    lesserG = cRibbon_K.LesserGreenFunction();
    
    % check the equivalence with the equlibrium lesser Green operator
    retardedG = cRibbon_K.FL_handles.Read('G');
    occupation_number = cRibbon_K.Fermi( Energy );
    disp(['Occupation number = ', num2str( occupation_number ), ', Energy: ', num2str(Energy), ', temperature = ', num2str(T)]);
    
    checkpoint = norm( lesserG - occupation_number*(retardedG'-retardedG ) );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
    
    %% test of class NTerminal_Keldysh, Transport_Interface_Keldysh
    disp( '***** test of class NTerminal_Keldysh, Transport_Interface_Keldysh *****' );
    
    % create input structures for three-terminal setup
    [Opt, param] = parseInput( 'Input_ThreeTerminal.xml');
    
    % chemical potentials in the leads
    bias_leads = [0.0 0.0 0.0];
    
    % chemical potential in the central device
    EF = 0.05;
    
    % creating function handle for the custom Hamiltonians
    hHamiltonians = @(varargin)ThreeTerminalHamiltonians(Opt, param, varargin{:});
    
    % set new temperature value
    T = 1000;  
    
    % Set the Energy value
    Energy = 0.04;
    
    % create an instance of class Lead_Keldysh
    cNTerminal_K = NTerminal_Keldysh('Opt', Opt, 'param', param, 'E', Energy, 'bias_leads', bias_leads, 'EF', EF, 'T', T, 'CustomHamiltoniansHandle', hHamiltonians);    
    
    % calculate the inverse of the retarded surface Green operator of the isolated scattering center
    cNTerminal_K.CalcFiniteGreensFunction('onlyGinv', true);
    
    % calculate the lesser Green operator
    lesserG = cNTerminal_K.LesserGreenFunction();
    
    % check the equivalence with the equlibrium lesser Green operator
    retardedG = cNTerminal_K.FL_handles.Read('G');
    occupation_number = cNTerminal_K.Fermi( Energy );
    disp(['Occupation number = ', num2str( occupation_number ), ', Energy: ', num2str(Energy), ', temperature = ', num2str(T)]);
    
    checkpoint = norm( lesserG - occupation_number*(retardedG'-retardedG ) );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
    
    
    
%% test of class NTerminal_Keldysh for biased leads
    disp( '***** test of class NTerminal_Keldysh, biased leads *****' );
    
    % create input structures for three-terminal setup
    [Opt, param] = parseInput( 'Input_ThreeTerminal.xml');
    
    % chemical potentials in the leads
    bias_leads = [0.0 0.0 0.1];
    
    % chemical potential in the central device
    EF = 0.05;
    
    % creating function handle for the custom Hamiltonians
    hHamiltonians = @(varargin)ThreeTerminalHamiltonians(Opt, param, varargin{:});
    
    % set new temperature value
    T = 10;  
    
    % Set the Energy value
    Energy = 0.04;
    
    % create an instance of class Lead_Keldysh
    cNTerminal_K = NTerminal_Keldysh('Opt', Opt, 'param', param, 'E', Energy, 'bias_leads', bias_leads, 'EF', EF, 'T', T, 'CustomHamiltoniansHandle', hHamiltonians);    
    
    % calculate the inverse of the retarded surface Green operator of the isolated scattering center
    cNTerminal_K.CalcFiniteGreensFunction('onlyGinv', true);
    
    % calculate the lesser Green operator
    [lesserG, junction_sites] = cNTerminal_K.LesserGreenFunction();
    
    % **** check the Dyson equation of Eq (6) of Eur. Phys. J. B 53, 537-549 (2006) ****
        
    % Calculating the retarded and lesser Green operator of the uncoupled system
    retardedg = zeros(size(lesserG));
    lesserg = zeros(size(lesserG));
        
        
    % Calculating the Green operator of the isolated scattering center
    cNTerminal_K.CalcFiniteGreensFunctionFromHamiltonian();
    gfin_central = cNTerminal_K.GetFiniteGreensFunction();
    retardedg(end-size(gfin_central,1)+1:end, end-size(gfin_central,2)+1:end) = gfin_central;
        
    % calculating the occupation number of the central region
    occupation_number = cNTerminal_K.Fermi( Energy );
    lesserg(end-size(gfin_central,1)+1:end, end-size(gfin_central,2)+1:end) = occupation_number*(gfin_central'-gfin_central);
        
    % Calculating the Green operator of the leads
    Leads = cNTerminal_K.FL_handles.Read( 'Leads' );
    offset = 0;
    for idx = 1:length(Leads)
        Leads{idx}.SurfaceGreenFunction( 'CalcInverse', false );            
        gsurf = Leads{idx}.Read('gsurf');            
        retardedg( offset+(1:size(gsurf,1)), offset+(1:size(gsurf,2)) ) = gsurf;
            
        Leads{idx}.LesserSurfaceGreenFunction();
        lesserg_tmp = Leads{idx}.Read('lessergsurf');
        lesserg( offset+(1:size(gsurf,1)), offset+(1:size(gsurf,2)) ) = lesserg_tmp;
            
        offset = offset + size(gsurf,1);
    end
        
        
    % obtaining the coupling Hamiltonian
    Hcoupling = zeros( size(lesserG) );
    offset = 0;
    offset_adj = 0;
    for idx = 1:length(Leads)
        cNTerminal_K.CreateInterface( idx );
        tmp = cNTerminal_K.Interface_Regions{idx}.Read('Kcoupling');
        Hcoupling( offset+(1:size(tmp,1)), end-size(tmp,2)+1:end ) = tmp;
        offset = offset + size(tmp,1);
            
        tmp = cNTerminal_K.Interface_Regions{idx}.Read('Kcouplingadj');
        Hcoupling( end-size(tmp,1)+1:end, offset_adj+(1:size(tmp,2)) ) = tmp;
        offset_adj = offset_adj + size(tmp,2);
            
            
    end
        
        
    lesserG2 = lesserg + retardedG*Hcoupling*lesserg + lesserG*Hcoupling*retardedg';
    checkpoint = norm( lesserG2 - lesserG );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
        
        
    % **** checkpoint 3: check the Dyson equation of Eq (13) of PRB 57 7366 ****
        
    % create a global matrix containing the retarded self-energies and obtaining the retarded self energies of the leads
    retarded_self_energy = zeros( size(retardedG) );
    Leads = cNTerminal_K.FL_handles.Read( 'Leads' );
    leadnum = length( Leads );
    for idx = 1:leadnum
        indexes = junction_sites.Leads{idx}.site_indexes;
        retarded_self_energy( indexes, indexes ) = Leads{idx}.SelfEnergy();
    end
        
    lesserG3 = (eye(size(lesserG)) + retardedG*Hcoupling)*lesserg*(eye(size(lesserG)) + Hcoupling'*retardedG');
    checkpoint = norm( lesserG3 - lesserG );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end    
        
    % **** checkpoint 4: checking the Dyson equation of the retarded Green operator
    retardedG_2 = retardedg + retardedg*Hcoupling*retardedG;
    checkpoint = norm( retardedG_2 - retardedG );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end     
    


end