Self_Consistent_Spectrum.m

Go to the documentation of this file.

%%    Calculates the self-consitent spectrum of a graphene quantum dot
%    Compare results to Fig.2 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 Self_Consistent_Spectrum.m
%> @brief Example to calculate the self-consitent spectrum of a graphene quantum dot
%> @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_Spectrum.jpg
%> @image latex Self_Consistent_Spectrum.jpg
%> @}
%
%> @brief Example to calculate the self-consitent spectrum of a graphene quantum dot
%> @param filenum The identification number of the filenema for the exported data (default is 1).
function Self_Consistent_Spectrum( filenum )

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

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

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

    % Strength of the self-consitent intercation in the Hubbard model
    U = -0.6*vargamma;
    % Zeeman energy for the initial iteration to bring the system out of the trivial solution
    Zeeman_energy = 0.5*vargamma; 
    % The Fermi energy
    EF = 0; % in eV
    % 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;
    
    % the calculated energy eigenvalues in the non-interacting and in the self-consistent cases
    filling_factors = [];
    eig_spinup = [];
    eig_no_interaction =  [];
    
    
    % Class for self-consistent iterations
    cSelfConsistent  = [];
    
    % creating output directories
    setOutputDir();   
    
    % Creates the class for self-consistent iterations    
    cSelfConsistent = SelfConsistent( Opt, 'magnetizaion_direction', []) ; 
    
    % create the an instance of class Ribbon describing the structure
    cRibbon = Ribbon('width', width, 'height', height, 'EF', EF, 'silent', silent, 'Opt', Opt, 'param', param, 'filenameOut', fullfile(outputdir, [outfilename, '.xml']) );
    
    % Calculates the self-consistent spectrum of graphene island
    CalcSelfConsistentSpectrum()
    
    % creates the plot
    PlotFunction()

    % export the calculeted data
    save( [outputdir, filesep, outfilename,'.mat']);  
    
    
%% CalcSelfConsistentSpectrum    
%> @brief Calculates the self-consistent spectrum of a graphene island
    function CalcSelfConsistentSpectrum()
        
        % creating an array of band filling factors for which the
        %calculation shall be performed
        filling_factor = floor(width*(2*height-1)/2)*2;
        filling_factors =  filling_factor-8:2:filling_factor+10;
        
        % preallocating arrays for the eienenergies
        eig_num = length(filling_factors);       
        eig_spinup = zeros(eig_num,1);
        eig_no_interaction =  zeros(eig_num,1);    
    
        % calculating the non-interacting eigenvalues
        
        % creating the Hamiltonian of the scattering region
        cRibbon.CreateScatter();
        CreateH = cRibbon.CreateH;
        ScatterPot(CreateH, 0); %removing the unecessary sites from the Hamiltonian
        Hscatter = CreateH.Read('Hscatter');
        eig_no_interaction = eig(Hscatter);
        eig_no_interaction = eig_no_interaction(filling_factors)/vargamma;
        
        % calculating the self-consistent eigenvalues
        
        % calculate the sel-consistent potential
        CalcSelfconsistentPotential( filling_factor )
        
        % creating the Hamiltonian of the scattering region
        cRibbon.CreateScatter();
        CreateH = cRibbon.CreateH;
        pot = ScatterPot(CreateH, 0);
        Hscatter = CreateH.Read('Hscatter');
        Hscatter_Self_consistent = Hscatter + diag(pot); 
            
        % calculate the self-consistent energy eigenvalues
        eigvals = eig(Hscatter_Self_consistent)/vargamma;         
        eig_spinup = eigvals(filling_factors);   
        
        
    end
    

    
%% CalcSelfconsistentPotential
%> @brief Calculates the self-consistent potential of a graphene island
    function CalcSelfconsistentPotential( filling_factor )     
        % reset the self-consistent class to forget previous iterations, but keep the last calculated density
        cSelfConsistent.Reset();
        
        % Do the self-consistent iterations with function handle pointing to a Stoner model
        cSelfConsistent.runSelfConsistentIterations( @()DensityCalc(filling_factor) );
   
        
    end
    
    
    
%% DensityCalc
%> @brief Calculates the density using the predefined graphene lattice framework (Ribbon class)
    function [density_data, junction_sites_spinup] = DensityCalc( filling_factor )
                
        % create Density class
        cDensity = Density(Opt, 'T', T, 'scatterPotential', @(CreateH, E)(ScatterPot(CreateH, E)), 'junction', cRibbon);
        
        % Calculate the density   
        [density_data, junction_sites_spinup] = cDensity.DensityCalcSmall(filling_factor); 
    
    end

%% ScatterPot
%> @brief Stoner model in the scattering region (and removing unnecessary sites from the scattering region) 
    function ret = ScatterPot(CreateH, E)
        
        
        % determine sites to be removed to reproduce the scattering region in  2014 Nanotechnology 25 465201
        CreateH.Write('kulso_szabfokok', []);
          
        coordinates = CreateH.Read('coordinates');
        y = coordinates.y;
        max_y = max(y);
            
        indexes2remove = ~(y<max_y);
        CreateH.RemoveSites( indexes2remove )        

        
        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(~coordinates.spinup) = -ret(~coordinates.spinup);
            return     
        end
        
        % Hubbard model
%         if isempty(cSelfConsistent.density_spinup) && spin==1
%             density_scatter_spindown =cSelfConsistent. density_spindown(junction_sites_spindown.Scatter.site_indexes); 
%              ret = ret + U*(density_scatter_spindown);       
%         elseif  isempty(cSelfConsistent.density_spindown) && spin==-1
%             density_scatter_spinup = cSelfConsistent.density_spinup(junction_sites_spinup.Scatter.site_indexes);
%              ret = ret + U*(density_scatter_spinup);     
%         end        
        
        % Stoner interaction model of Nanotechnology 25 (2014) 465201
        coordinates = CreateH.Read('coordinates');
        ret = zeros(size(coordinates.x));
        
        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



%% PlotFunction
%> @brief Creates the plot
    function PlotFunction()
        
        % creating figure in units of pixels
        figure1 = figure( 'Units', 'Pixels', 'Visible', 'off');
        
        % font size on the figure will be 16 points
        fontsize = 16;         
        
        % creating axes of the plot        
        axes1 = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,...
                'Box', 'on',...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        % Create xlabel
        xlabel('Band filling','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);

        % Create ylabel
        ylabel('E [\gamma]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes1);
        
        % reserve a variable for legend titles
        legend_labels = [];
        
        % plot the non-interacting and self-consistent data
        number_of_states = width*(2*height-1)*2;
        if ~isempty( eig_no_interaction )
            plot(filling_factors/number_of_states, eig_no_interaction, 'bo', 'Parent', axes1)
               legend_labels{end+1} = 'U = 0';
        end
        
        if ~isempty( eig_spinup )
            plot(filling_factors/number_of_states, eig_spinup, 'rx', 'Parent', axes1)
               legend_labels{end+1} = ['U = ', num2str(U/vargamma), '\gamma spin up'];
        end
                
        % set the legends
        fig_legend = legend(axes1, legend_labels);
        set(fig_legend, 'FontSize', fontsize*0.8, 'FontName','Times New Roman', 'Box', 'off', 'Location', 'SouthEast')
        
        % 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
%> @brief Sets output directory.
    function setOutputDir()
        resultsdir = 'results';
        mkdir(fullfile(directory, resultsdir) );
        outputdir = resultsdir;        
    end

        
        
    
end