Coulomb_Diamonds.m

Go to the documentation of this file.

%    Transport calculation through a graphene QD including the charging energy - based on EQuUs v4.8
%    Compare results to Fig.3 in 2014 Nanotechnology 25 465201.
%    Copyright (C) 2017 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 Coulomb_Diamonds.m
%> @brief Transport calculation through a graphene QD including the charging energy utilizing the diagarmmatic technique of ......
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @reference
%> Herbert Scholler and Gerd Schon Phys RevB <B>50</B> 18436 (1994)
%> @Available
%> EQuUs v4.9 or later
%> @expectedresult
%> @image html Coulomb_Diamonds.jpg
%> @image latex Coulomb_Diamonds.jpg
%> @}
%
%> @brief Transport calculation through a graphene QD including the charging energy utilizing the diagarmmatic technique of ......
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function Coulomb_Diamonds( filenum )

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

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

    % filename containing the input XML    
    inputXML = 'Input_Coulomb.xml';
    % Parsing the input file and creating data structures
    [Opt, param] = parseInput( fullfile( directory, inputXML ) ); 
    
    % The hopping parameter obtained from the input file
    vargamma = param.scatter.vargamma;

    % The charging energy
    Ec = 0.1*vargamma;    
    % The maximal bias on the system
    bias_max = 4*Ec;
    % The maximal bias on the system
    gateVoltage_max = 5*Ec;    
    % number of energy points to be used for the calculations in the Quantum Dot class
     Enum = 800;    
    % The Boltzman constant
    k_B = 8.6173324e-5;    
    % temperature in K
    T = 0.05*Ec/k_B;
    % width of the junction    
    width = 6;
    % length of the junction
    height = 8;
    % filename containing the output XML
    outfilename = [fncname,'_',num2str( filenum )];
    % the output folder
    outputdir   = [];
    % false to supress output messages
    silent = false;    
    % The differential conductance as a function of the magnetic field
    diff_conductance = [];
    
     % array containing the bias voltage values
    bias_vec = [];
     % array containing the gate voltage values
    gateVoltage_vec = [];
     % Ribbon class describing the junction
    cRibbon = [];
     % Two dimensional array containg the calculated current values
     current_surf = [];
    
    % doping levels of the QD in the single-electron picture (without the charging energy)
    %doping_levels = [ -0.981 -0.963 -0.276 -0.246 -0.009179 0  0.009179 0.036 0.2754 0.309 0.981 1.00]; %gamma_sc = 0.3    
    doping_levels = [ -0.993 -0.984 -0.981 -0.951 -0.279 -0.231 -0.009179 0  0.009179 0.057 0.279 0.342 0.981 0.984 0.993 1.026]; %gamma_sc = 0.5   
    

    
    % creating output directories
    setOutputDir();    
        
    % Calculates the non-interacting transport through a graphene island
    CalcTransport()
    
    % creates the plot
    PlotFunction()
    
    % exporting the calculated data
    save( fullfile(outputdir, [outfilename, '.mat']));
    
%% CalcTransport
%> @brief Calculates the transport through a quantum dot
    function CalcTransport()
        
        display('****** Calculating non-interacting transport through a Quantum dot ******')
        display(' ')   
        
        % creating the Ribbon interface describing the junction
        cRibbon = Ribbon('width', width, 'height', height, 'EF', 1e-6, 'silent', silent, 'Opt', Opt, 'param', param, 'filenameOut', fullfile(outputdir, [outfilename, '.xml']), ...
            'leadmodel', @Square_leads, 'interfacemodel', @InterfaceModel);        
        
            
        
        % creating array containing the bias values
        bias_num = 200;
        bias_min = -bias_max;        
        bias_vec = bias_min:(bias_max-bias_min)/(bias_num+1):bias_max;  
        bias_vec = sort( [bias_vec, 0], 'ascend');
        
        % creating array containing the gate voltage values
        gateVoltage_num = 200;
        gateVoltage_min = -gateVoltage_max;
        gateVoltage_vec = gateVoltage_min:(gateVoltage_max-gateVoltage_min)/(gateVoltage_num+1):gateVoltage_max + Ec;
             
        
        % creating the QuantumDot class to calculate the current
        cQD = QuantumDot( Opt, 'Ec', Ec, 'T', T, 'bias_max', bias_max, 'junction', cRibbon, 'scatterPotential', @(CreateH,E)ScatterPotQD(CreateH, E),...
            'gfininvfromHamiltonian', true, 'Edb', Enum);
        
        
        % preallocating arrays for the calculations  
        current_surf = zeros( length(gateVoltage_vec), length(bias_vec) );
        diff_conductance = zeros( length(gateVoltage_vec), length(bias_vec)-1 );

        
        % starting parallel pool
        poolmanager = Parallel( Opt );
        poolmanager.openPool();  
        
        % Determine the DOS for the calculations
        %cQD.DetermineDensityOfStates( -1.5, 1.5, 1000);
        
        for jdx = 1:length(gateVoltage_vec)
            
            gateVoltage = gateVoltage_vec(jdx);
            disp(['Calculating for gate voltage = ', num2str(gateVoltage)]); 
            
            % Determine chamical potentials of the individual filling factors
            cQD.SetChemicalPotentials( doping_levels, gateVoltage ); 

            % calculating the tunneling rates for the given gate voltage
            cQD.CalcTunnelingRates( gateVoltage );

        
            for idx = 1:length(bias_vec)
        
                % the current bias
                bias = bias_vec(idx);  

                % determine the transition rates 
                cQD.DetermineSelfEnergy( bias, gateVoltage );
            
                % determine the steady state occupation probabilities
                OccupationProbabilities = cQD.DetermineOccupationProbabilities();


                    if abs(bias - min(abs(bias_vec))) < 1e-4
                         display('Occupation probabilities in equilibrium')
                         disp(transpose(OccupationProbabilities));
                    end
            
                % Determine the new current operator for the new bias value at lead 1
                cQD.DetermineCurrentOperator( bias, gateVoltage, 1 );        
            
                % calculate the steady state current
                    current_tmp = cQD.DetermineCurrent();   
                current_surf(jdx,idx) = current_tmp;

                    % Check the equivalence of the current in the two leads
                    if abs(bias - max(bias_vec)) < 1e-4

                         % Determine the new current operator for the new bias value at lead 2
                     cQD.DetermineCurrentOperator( bias, gateVoltage, 2 );        
            
                    % calculate the steady state current
                     current_tmp2 = cQD.DetermineCurrent();  

                         disp(['The difference between the current of the two leads: ', num2str( current_tmp+current_tmp2)])

                    end

            end                
           
            % calculating differential conductance
            diff_conductance(jdx, :) = diff(current_surf(jdx,:))./diff(bias_vec);
            
            % creates the plot
            PlotFunction()

          % export the calculetd data
          save( [outputdir,'/',outfilename,'.mat']);  
            
        end       
              
        % closing the parallel pool
        poolmanager.closePool();                
           
            
                
    end    



%% ScatterPotQD
%> @brief Removing unecessary sites from the scattering region
    function ret = ScatterPotQD(CreateH, E )
        
        coordinates = CreateH.Read('coordinates');
        y = coordinates.y;
        max_y = max(y);
        
        % determine sites to be removed to reproduce the scattering region in  2014 Nanotechnology 25 465201
        indexes2remove = ~(y<max_y);
        
        CreateH.Write('kulso_szabfokok', []);
        
        y( indexes2remove ) = max_y/2;
        max_y = max(y);
        min_y = min(y);
        
        % creating new array for the edge-sites
        edge_sites_logical = ~(y<max_y & y>min_y);
        edge_sites = 1:length(y);
        edge_sites = edge_sites(edge_sites_logical);
        CreateH.Write('kulso_szabfokok', edge_sites);
        

        % removing sites from the scattering region
        CreateH.RemoveSites( indexes2remove )        
        
        
        % store the modified Hamiltonian
        Hscatter = CreateH.Read('Hscatter'); 
        
        ret = sparse([], [], [], size(Hscatter,1), size(Hscatter,2));
        
    end

%% Square_leads
%> @brief Function producing custom lead Hamiltonians defined on a square lattice. The function has equivalent inputs and return values as #Transport_Interface.SurfaceGreenFunctionCalculator and E is standing for the energy.
    function Lead_ret = Square_leads( idx, E, varargin )
        p_inner = inputParser;
        
        p_inner.addParameter('createCore', 0);
        p_inner.addParameter('Just_Create_Hamiltonians', 0);
        p_inner.addParameter('shiftLead', 0);
        p_inner.addParameter('coordinates_shift', 0);
        p_inner.addParameter('transversepotential', []);
        p_inner.addParameter('Lead',[]);
        p_inner.addParameter('gauge_field', [] );% gauge field for performing gauge transformation
        p_inner.addParameter('SelfEnergy',false); % set true to calculate the self energy of the semi-infinite lead
        p_inner.addParameter('SurfaceGreensFunction', true );% set true to calculate the surface Greens function of the semi-infinite lead
        p_inner.addParameter('q', [] ) %transverse momentum
        
        p_inner.parse(varargin{:});
        
        createCore            = p_inner.Results.createCore;
        Just_Create_Hamiltonians = p_inner.Results.Just_Create_Hamiltonians;
        shiftLead             = p_inner.Results.shiftLead;
        coordinates_shift     = p_inner.Results.coordinates_shift;
        transversepotential   = p_inner.Results.transversepotential;
        Lead                  = [];%p_inner.Results.Surface_tmp;
        gauge_field           = p_inner.Results.gauge_field; % gauge field for performing gauge transformation
        SelfEnergy            = p_inner.Results.SelfEnergy;
        SurfaceGreensFunction = p_inner.Results.SurfaceGreensFunction;
        q                     = p_inner.Results.q;
        
        Opt2 = Opt;
        Opt2.Lattice_Type = 'Square';
        Opt2.Silent = true;
        
        param2 = param;
        param2.Leads{idx} = param_Square_Lead();
        param2.Leads{idx}.epsilon = param.Leads{idx}.epsilon;
        param2.Leads{idx}.vargamma = param.Leads{idx}.vargamma;
        param2.Leads{idx}.M = width;
        
        
        FL_handles = Transport_Interface(E, Opt2, param2);
        
        Lead_ret = FL_handles.SurfaceGreenFunctionCalculator(idx, 'createCore', createCore, ...
                            'Just_Create_Hamiltonians', Just_Create_Hamiltonians, ...
                            'shiftLead', shiftLead, ...
                            'coordinates_shift', coordinates_shift, ...
                            'transversepotential', transversepotential, ...
                            'gauge_field', gauge_field, ...
                            'SelfEnergy', SelfEnergy, ...
                            'SurfaceGreensFunction', SurfaceGreensFunction, ...
                            'q', q);
        
    end


%% InterfaceModel
%> @brief Method to adjust the Hamiltonians of the interface regions
    function InterfaceModel( Interface_Region )
        
        Kcoupling = Interface_Region.Read('Kcoupling');
        Kcouplingadj = Interface_Region.Read('Kcouplingadj');
            
            
        [rows, cols] = find( Kcoupling );
        rows = unique(rows);          
            
        Kcoupling = Kcoupling(rows, :);
        Kcouplingadj = Kcouplingadj(:, rows);
            
        Interface_Region.Write('Kcoupling', Kcoupling);
        Interface_Region.Write('Kcouplingadj', Kcouplingadj); 

        
    end

%% PlotFunction
%> @brief Creates the plot
    function PlotFunction()
            
        % creating figure in units of pixels
        figure1 = figure('Units', 'Pixels', 'Visible', 'off', 'Colormap', ...
    [0 0 0;0.0300625283271074 0.0187878794968128 0.0119648873806;0.0601250566542149 0.0375757589936256 0.0239297747612;0.0901875868439674 0.0563636384904385 0.0358946621417999;0.12025011330843 0.0751515179872513 0.0478595495223999;0.150312647223473 0.0939393937587738 0.0598244369029999;0.180375173687935 0.112727276980877 0.0717893242835999;0.210437700152397 0.131515145301819 0.0837542116641998;0.240500226616859 0.150303035974503 0.0957190990447998;0.270562767982483 0.169090911746025 0.1076839864254;0.300625294446945 0.187878787517548 0.119648873806;0.330687820911407 0.20666666328907 0.1316137611866;0.36075034737587 0.225454553961754 0.1435786485672;0.390812873840332 0.244242429733276 0.1555435359478;0.420875400304794 0.263030290603638 0.1675084233284;0.450937926769257 0.281818181276321 0.179473310709;0.481000453233719 0.300606071949005 0.1914381980896;0.511062979698181 0.319393932819366 0.2034030854702;0.541125535964966 0.33818182349205 0.2153679728508;0.571188032627106 0.356969714164734 0.2273328602314;0.60125058889389 0.375757575035095 0.239297747612;0.63131308555603 0.394545465707779 0.251262634992599;0.661375641822815 0.41333332657814 0.263227522373199;0.6914381980896 0.432121217250824 0.275192409753799;0.72150069475174 0.450909107923508 0.287157297134399;0.751563251018524 0.469696968793869 0.299122184514999;0.781625747680664 0.488484859466553 0.311087071895599;0.811688303947449 0.507272720336914 0.323051959276199;0.841750800609589 0.526060581207275 0.335016846656799;0.871813356876373 0.544848501682281 0.346981734037399;0.901875853538513 0.563636362552643 0.358946621417999;0.931938409805298 0.582424223423004 0.370911508798599;0.962000906467438 0.60121214389801 0.382876396179199;0.992063462734222 0.620000004768372 0.394841283559799;0.992327988147736 0.625373363494873 0.398263245820999;0.992592573165894 0.63074666261673 0.401685208082199;0.992857098579407 0.636120021343231 0.405107140541077;0.993121683597565 0.641493320465088 0.408529102802277;0.993386209011078 0.646866679191589 0.411951065063477;0.993650794029236 0.652239978313446 0.415373027324677;0.993915319442749 0.657613337039948 0.418794989585876;0.994179844856262 0.662986695766449 0.422216951847076;0.99444442987442 0.668359994888306 0.425638914108276;0.994708955287933 0.673733353614807 0.429060846567154;0.994973540306091 0.679106652736664 0.432482808828354;0.995238065719604 0.684480011463165 0.435904771089554;0.995502650737762 0.689853310585022 0.439326733350754;0.995767176151276 0.695226669311523 0.442748695611954;0.996031761169434 0.700600028038025 0.446170628070831;0.996296286582947 0.705973327159882 0.449592590332031;0.99656081199646 0.711346685886383 0.453014552593231;0.996825397014618 0.71671998500824 0.456436514854431;0.997089922428131 0.722093343734741 0.459858477115631;0.997354507446289 0.727466642856598 0.463280439376831;0.997619032859802 0.732840001583099 0.466702401638031;0.99788361787796 0.738213300704956 0.470124334096909;0.998148143291473 0.743586659431458 0.473546296358109;0.998412668704987 0.748960018157959 0.476968258619308;0.998677253723145 0.754333317279816 0.480390220880508;0.998941779136658 0.759706676006317 0.483812183141708;0.999206364154816 0.765079975128174 0.487234115600586;0.999470889568329 0.770453333854675 0.490656077861786;0.999735474586487 0.775826632976532 0.494078040122986;1 0.781199991703033 0.497500002384186]);

        % 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('Gate voltage [E_c]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);

        % Create ylabel
        ylabel('Bias [E_c]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);
                
        % plot the data
        if ~isempty( diff_conductance )
            contourf( gateVoltage_vec/Ec, (bias_vec(1:end-1)+bias_vec(2:end))/2/Ec,  transpose(diff_conductance), 'LineStyle','none','LineColor',[0 0 0],...
    'LevelList',[-0.000275211478671962 0.00015807932351996 0.00059137012571188 0.0010246609279038 0.00145795173009572 0.00189124253228764 0.00232453333447956 0.00275782413667149 0.00319111493886341 0.00362440574105533 0.00405769654324725 0.00449098734543917 0.00492427814763109 0.00535756894982301 0.00579085975201493 0.00622415055420685 0.00665744135639878 0.0070907321585907 0.00752402296078262 0.00795731376297454 0.00839060456516646 0.00882389536735838 0.0092571861695503 0.00969047697174222 0.0101237677739341 0.0105570585761261 0.010990349378318 0.0114236401805099 0.0118569309827018 0.0122902217848937 0.0127235125870857 0.0131568033892776 0.0135900941914695 0.0140233849936614 0.0144566757958534 0.0148899665980453 0.0153232574002372 0.0157565482024291 0.016189839004621 0.016623129806813 0.0170564206090049 0.0174897114111968 0.0179230022133887 0.0183562930155806 0.0187895838177726 0.0192228746199645 0.0196561654221564 0.0200894562243483 0.0205227470265402 0.0209560378287322 0.0213893286309241 0.021822619433116 0.0222559102353079 0.0226892010374999 0.0231224918396918 0.0235557826418837 0.0239890734440756 0.0244223642462675 0.0248556550484595 0.0252889458506514 0.0257222366528433 0.0261555274550352 0.0265888182572271 0.0270221090594191 0.027455399861611 0.0278886906638029 0.0283219814659948 0.0287552722681867 0.0291885630703787 0.0296218538725706 0.0300551446747625 0.0304884354769544 0.0309217262791464 0.0313550170813383 0.0317883078835302 0.0322215986857221 0.032654889487914 0.033088180290106 0.0335214710922979 0.0339547618944898 0.0343880526966817 0.0348213434988736 0.0352546343010656 0.0356879251032575 0.0361212159054494 0.0365545067076413 0.0369877975098332 0.0374210883120252 0.0378543791142171 0.038287669916409 0.0387209607186009 0.0391542515207929 0.0395875423229848 0.0400208331251767 0.0404541239273686 0.0408874147295605 0.0413207055317525 0.0417539963339444 0.0421872871361363 0.0426205779383282 0.0430538687405201 0.0434871595427121 0.043920450344904 0.0443537411470959 0.0447870319492878 0.0452203227514798 0.0456536135536717 0.0460869043558636 0.0465201951580555 0.0469534859602474 0.0473867767624394 0.0478200675646313 0.0482533583668232 0.0486866491690151 0.049119939971207 0.049553230773399 0.0499865215755909 0.0504198123777828 0.0508531031799747 0.0512863939821666 0.0517196847843586 0.0521529755865505 0.0525862663887424 0.0530195571909343 0.0534528479931263 0.0538861387953182 0.0543194295975101 0.054752720399702 0.0551860112018939 0.0556193020040859 0.0560525928062778 0.0564858836084697 0.0569191744106616 0.0573524652128535 0.0577857560150455 0.0582190468172374 0.0586523376194293 0.0590856284216212 0.0595189192238131 0.0599522100260051 0.060385500828197 0.0608187916303889 0.0612520824325808 0.0616853732347728 0.0621186640369647 0.0625519548391566 0.0629852456413485 0.0634185364435404 0.0638518272457324 0.0642851180479243 0.0647184088501162 0.0651516996523081 0.0655849904545 0.066018281256692 0.0664515720588839 0.0668848628610758 0.0673181536632677 0.0677514444654596 0.0681847352676516 0.0686180260698435 0.0690513168720354 0.0694846076742273 0.0699178984764192 0.0703511892786112 0.0707844800808031 0.071217770882995 0.0716510616851869 0.0720843524873788 0.0725176432895708 0.0729509340917627 0.0733842248939546 0.0738175156961465 0.0742508064983385 0.0746840973005304 0.0751173881027223 0.0755506789049142 0.0759839697071061 0.0764172605092981 0.07685055131149 0.0772838421136819 0.0777171329158738 0.0781504237180657 0.0785837145202577 0.0790170053224496 0.0794502961246415 0.0798835869268334 0.0803168777290253 0.0807501685312173 0.0811834593334092 0.0816167501356011 0.082050040937793 0.082483331739985 0.0829166225421769 0.0833499133443688 0.0837832041465607 0.0842164949487526 0.0846497857509446 0.0850830765531365 0.0855163673553284 0.0859496581575203 0.0863829489597123],...
    'Fill','on', 'Parent', axes1);
        end
        
        % 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
              
        Position_figure = get(figure1, 'Position');
        Position_figure(1) = 0;
        Position_figure(2) = 0;
        set(axes1, 'OuterPosition', Position_figure);

        % 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
         
    

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

        
        
    
end