Self_Consistent_Transport.m

Go to the documentation of this file.

%%    Transport calculation through a graphene QD using self-consistent potential - based on EQuUs v4.8
%    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 Self_Consistent_Transport.m
%> @brief Example to calculate transport through a graphene quantum dot using self-consistent potential. Compare results to Fig.3 in 2014 Nanotechnology 25 465201. In these calculations we use a leads defined on a square lattice and do the self-consistent potential iterations only for \f$E_F =0\f$. 
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @reference
%> S. Krompiewski, Nanotechnology 25 465201 (2014).
%> @Available
%> EQuUs v4.9 or later
%> @expectedresult
%> @image html Self_Consistent_Transport.jpg
%> @image latex Self_Consistent_Transport.jpg
%> @}
%
%> @brief Example to calculate transport through a graphene quantum dot using self-consistent potential. Compare results to Fig.3 in 2014 Nanotechnology 25 465201. In these calculations we use a leads defined on a square lattice and do the self-consistent potential iterations only for \f$E_F =0\f$. 
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function Self_Consistent_Transport( filenum )

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

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

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

    % Strength of the self-consitent intercation in the Stoner model
    U = -0.6*vargamma;    
    % Zeeman energy for the initial iteration to bring the system out of the trivial nonmagnetic solution
    Zeeman_energy = 0.1*vargamma;   
    % The Fermi energy
    EF = 1e-6; % in eV    
    % number of energy point on the contour
    Edb = 51;    
    % temperature in K
    T = 0;
    % width of the junction    
    width = 22;
    % length of the junction
    height = 9;
    % filename containing the output XML
    outfilename = [fncname,'_',num2str( filenum )];
    % the output folder
    outputdir   = [];
    % false to supress output messages
    silent = false;
    % Using self energy in the transport calculations
    useselfEnergy = true;
    
    
    % noninteracting transport results
    Evec_non_interacting = [];
    Conductance_non_interacting = [];
    
    % self_consistent transport results
    Conductance_self_consistent = [];
    Evec_self_consistent = [];
    
    % Class for self_consistent iterations
    cSelfConsistent  = [];
    
    % creating output directories
    setOutputDir();    
        
    % Calculates the non-interacting transport through a graphene island
    CalcTransport()
    
    % Calculates the self-consistent transport through a graphene island
    CalcSelfConsistentTransport()    
    
    % creates the plot
    CreatePlot()
    
    % export the calculetd data
    save( [outputdir,'/',outfilename,'.mat']);  
    
%% CalcTransport
%> @brief calculates the transport with self-consistent potential
    function CalcSelfConsistentTransport()
        
        display('****** Calculating self-consistent transport ******')
        display(' ')
        
        % Creates the class for self-consistent iterations for magnetization up
        cSelfConsistent = SelfConsistent( Opt, 'magnetizaion_direction', [], 'tolerance', 2e-3, 'newIterationWeight', 0.3, 'outFileName', 'test.mat', 'MaxIter', 500, 'plotIterations', false, 'resume', false ) ; 
        
        Enum = 100;
        Emin = -0.3;
        Emax = 0.3;
        Evec_self_consistent = Emin:(Emax-Emin)/Enum:Emax;
        Conductance_self_consistent = zeros(size(Evec_self_consistent));
        
        % calcualting the self-consistent potential
        
        % opening the parallel pool
        Opt.workers = 2;
        poolmanager = Parallel( Opt );
        poolmanager.openPool();
        
        % running the self-consistent iterations
        CalcSelfconsistentPotential();
        
        % closing the parallel pool
        poolmanager.closePool();  
        
        for idx = 1:length(Evec_self_consistent)
            
            Conductance_self_consistent(idx) = Transport(Evec_self_consistent(idx), true);
            
            display(['Energy = ', num2str(Evec_self_consistent(idx)), ', conductance = ', num2str(Conductance_self_consistent(idx))]);           
    
            
        end
        
        
    end    
    
    
%% CalcTransport
%> @brief calculates the non-interacting transport
    function CalcTransport()
        
        display('****** Calculating transport in the non-interacting case ******')
        display(' ')
        
        % oreallocating arrays for the results
        Enum = 100;
        Emin = -0.3;
        Emax = 0.3;
        Evec_non_interacting = Emin:(Emax-Emin)/Enum:Emax;
        Conductance_non_interacting = zeros(size(Evec_self_consistent));
        
        for idx = 1:length(Evec_non_interacting)
            Conductance_non_interacting(idx) = Transport(Evec_non_interacting(idx), false);
            display(['Energy = ', num2str(Evec_non_interacting(idx)), ', conductance = ', num2str(Conductance_non_interacting(idx))])  
        end
        
        % creates the plot
        CreatePlot()        
        
        
    end


%% Transport
%> @brief Calculates the conductance for a given energy for both self-consistent and non-interacting cases.
%> @param Energy The enrgy value
%> @param self_consistent Logical value. Set true to use the self-consistent potential, or false to do the non-interacting calculations.
%> @return Returns with the calculated conductance in units of \f$e^2/\hbar \f$.
    function ret = Transport( Energy, self_consistent )
        
        
        
        % creating the Ribbon interface for the SNS structure
        cRibbon = Ribbon('width', width, 'height', height, 'EF', EF, 'silent', silent, 'Opt', Opt, 'param', param, 'filenameOut', fullfile(outputdir, [outfilename, '.xml']), ...
            'leadmodel', @Square_leads, 'interfacemodel', @InterfaceModel);
        
        % setting the energy value
        cRibbon.setEnergy( Energy )
        
        % calculating the Green function of the scattering reigon
        cRibbon.CalcFiniteGreensFunctionFromHamiltonian('PotInScatter', @(CreateH,E)Pot_Self_Consistent(CreateH, E, self_consistent), 'onlyGinv', true); 
        
        % calculating the leads
        coordinates_shift = [-2, height+1];
        recalculateLeads = [1 2];
          cRibbon.FL_handles.LeadCalc('coordinates_shift', coordinates_shift, 'SelfEnergy', useselfEnergy, 'SurfaceGreensFunction', ~useselfEnergy, 'leads', recalculateLeads, 'leadmodel', @Square_leads); 
        
        % creating interface regions
        for idx = 1:length(cRibbon.Interface_Regions)
            cRibbon.CreateInterface( idx );  

        end
        
        % setting the function handle for the Dyson Equation
        Dysonfunc = @()(cRibbon.CustomDysonFunc( 'recalculateSurface', [], 'constant_channels', true, 'decimate', true, 'SelfEnergy', useselfEnergy ));
        cRibbon.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
        
        % calculating the conductance 
        ret = cRibbon.FL_handles.Conductance2();
        
        
    end


%% CalcSelfconsistentPotential
%> @brief Calculates the self-consistent potential
    function CalcSelfconsistentPotential()     
        % reset the self-consistent class to forget previous iterations, but keep the last calculated density
        cSelfConsistent.Reset();             
        
        % Do the self-consistent iterations
        cSelfConsistent.runSelfConsistentIterations( @DensityCalc );
    end
    
    
    
%% DensityCalc
%> @brief Calculates the density using the predefined graphene lattice framework (Ribbon class)
%> @return [1] The calculated density data
%> @return [2] An instance of structure #junction_sites describing the sites in the calculated density
    function [density_data, junction_sites] = DensityCalc()        
       
        % creating the Ribbon interface for the SNS structure
        cRibbon = Ribbon('width', width, 'height', height, 'EF', EF, 'silent', silent, 'Opt', Opt, 'param', param, 'filenameOut', fullfile(outputdir, [outfilename, '.xml']), ...
            'leadmodel', @Square_leads, 'interfacemodel', @InterfaceModel);
        
        % create Density class
        cDensity = Density(Opt, 'T', T, 'scatterPotential', @(CreateH, E)(Pot_Self_Consistent(CreateH, E, true)), 'junction', cRibbon, 'useSelfEnergy', useselfEnergy);
        
        % Calculate the density   
        [density_data, junction_sites] = cDensity.DensityCalc( 'Edb', Edb, 'Emin', -11.0286 ); 
    
    end


%% Pot_Self_Consistent
%> @brief Create a potential from a self-consistent data, and removing unecessary sites from the scattering region
%> @param CreateH An instance of class #CreateHamiltonians storing the Hamiltonian.
%> @param E The energy value
%> @param SelfConsistent Logical value. Set true to calculate the self-consistent potential, or false otherwise
    function ret = Pot_Self_Consistent(CreateH, E, SelfConsistent)
        
        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 kulso_szabfokok
        kulso_szabfokok_logical = ~(y<max_y & y>min_y);
        kulso_szabfokok = 1:length(y);
        kulso_szabfokok = kulso_szabfokok(kulso_szabfokok_logical);
        CreateH.Write('kulso_szabfokok', kulso_szabfokok);
        

        % removing sites from the scattering region
        CreateH.RemoveSites( indexes2remove )        

        coordinates = CreateH.Read('coordinates');
        ret = zeros(size(coordinates.x));
        
        if ~SelfConsistent
            ret = zeros(size(coordinates.x));
            return
        end  
        
        
        % adding  the self-consistent interaction
        coordinates = CreateH.Read('coordinates');
        ret = zeros(size(coordinates.x));
        
        if isempty( cSelfConsistent.density)
            coordinates = CreateH.Read('coordinates');
            
            % adding Zeeman term at the edges the first iteration to bring the
            %system put of the trivial solution
            x = coordinates.x;
            min_x = min(x);
            max_x = max(x);
            lattice_constant = norm(coordinates.a);
            %ret = (exp(-(x-min_x).^2/(2*lattice_constant^2)) -  exp(-(x-max_x).^2/(2*lattice_constant^2)) )*Zeeman_energy;
            ret = -ones(size(x))*Zeeman_energy;
            ret(~coordinates.spinup) = -ret(~coordinates.spinup);
            return     
        end
        
        % Stoner interaction model of Nanotechnology 25 (2014) 465201
        scatter_density = cSelfConsistent.density(cSelfConsistent.junction_sites.Scatter.site_indexes);
        
        density_scatter_spinup = scatter_density(cSelfConsistent.junction_sites.Scatter.coordinates.spinup);
        density_scatter_spindown = scatter_density(~cSelfConsistent.junction_sites.Scatter.coordinates.spinup);
        
        ret(coordinates.spinup) = 0.5*U*(density_scatter_spinup-density_scatter_spindown);
        ret(~coordinates.spinup) = -0.5*U*(density_scatter_spinup-density_scatter_spindown);        
        
        
    end


%% Square lead
%> @brief Create square leads Hamiltonians with the same inputs and output as #Transport_Interface.SurfaceGreenFunctionCalculator.
    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
%> @param Interface_Region An instance of class #Interface_Region.
    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

%% CreatePlot
% creates the plot
    function CreatePlot()
        
        % creating figure in units of pixels
        figure1 = figure( 'Units', 'Pixels', 'Visible', 'off');
        
        % font size on the figure will be 16 points
        fontsize = 16;  
        
        % setting the limits of the axis
        x_lim = [-0.3 0.3];
        
        % creating axes of the plot        
        axes1 = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,...     
                'xlim', x_lim,...  
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        % Create xlabel
        xlabel('E [eV]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);

        % Create ylabel
        ylabel('G [e^2/h]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);
        
        % reserve a variable for legend titles
        legend_labels = [];
        
         % plot the non-interacting and self-consistent data
        if ~isempty( Conductance_non_interacting )
            plot(Evec_non_interacting, real(Conductance_non_interacting), 'b', 'LineWidth', 2, 'Parent', axes1)
               legend_labels{end+1} = 'U = 0';
        end
        
        if ~isempty( Conductance_self_consistent )
            plot(Evec_self_consistent, real(Conductance_self_consistent), 'rx', 'Parent', axes1)
            plot(Evec_self_consistent, real(Conductance_self_consistent), 'r', 'LineWidth', 1.5, 'Parent', axes1)
               legend_labels{end+1} = ['U = ', num2str(U/vargamma), '\gamma'];
        end
        
        
               
        % plot the non-interacting and self-consistent dat
        fig_legend = legend(axes1, legend_labels);
        set(fig_legend, 'FontSize', fontsize*0.8, 'FontName','Times New Roman', 'Box', 'off', 'Location', 'NorthEast')
    
        % 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
        
        set(axes1, 'Position', get(axes1, 'OuterPosition') - get(axes1, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1]);        
        Position_figure = get(axes1, 'OuterPosition');
        set(figure1, 'Position', Position_figure);
    

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

%% setOutputDir
% sets the output directory
    function setOutputDir()
        resultsdir = 'results';
        mkdir(fullfile(directory, resultsdir) );
        outputdir = resultsdir;        
    end

        
        
    
end