TMDC_Monolayer_Spectrum.m

Go to the documentation of this file.

%   Spectral density function in superconductor-graphene-superconductor junctions.
%    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 Examples
%> @{
%> @file TMDC_Monolayer_Spectrum.m
%> @brief Example to calculate the spectra of spinless monolayer \f$ M_oS_2 \f$ TMDC with nearest neighbor (EQuUs) and second nearest neighbor (<a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a>) model.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%> @Available
%> EQuUs v5.0 or later
%> @expectedresult
%> @image html TMDC_Monolayer_Spectrum.jpg
%> @image latex TMDC_Monolayer_Spectrum.jpg
%> @}
%
%> @brief Example to calculate the spectra of spinless monolayer \f$ M_oS_2 \f$ TMDC with nearest neighbor (EQuUs) and second nearest neighbor (<a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a>) model.
%> Test example for class TMDC_Monolayer_Lead_Hamiltonians.
%> @param filenum The identification number of the filenema for the exported data (default is 1).
%%
function TMDC_Monolayer_Spectrum( filenum )

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

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

    % filename containing the output XML
    outfilename = [fncname, '_',num2str( filenum )];
    % the output folder
    outputdir   = [];   
    
    % setting the output directory
     setOutputDir()

    % Parsing the input file and creating data structures
    [Opt, param] = parseInput( 'TMDC_Input.xml');

    % create an instance of class CreateLeadHamiltonians to create Hamiltonians
    HLead=CreateLeadHamiltonians(Opt, param, 'Hanyadik_Lead', 1);
    HLead.CreateHamiltonians();

    % getting lattice vectors of the triangle lattice
    cCoordinates = HLead.Read('coordinates');
    a1 = cCoordinates.a*cCoordinates.LatticeConstant;
    a2 = cCoordinates.b*cCoordinates.LatticeConstant;

    % nearest neighbour hopping vectors according to Table III in PRB 92 205108
    delta1 = a1;
    delta2 = a1 + a2;
    delta3 = a2;

    delta4 = -(2*a1+a2)/3;
    delta5 = (a1+2*a2)/3;
    delta6 = (a1-a2)/3;

    delta7 = -2*(a1+2*a2)/3;
    delta8 = 2*(2*a1+a2)/3;
    delta9 = 2*(a2-a1)/3;

    % reciprocal lattice vectors
    b1 = 2*pi/cCoordinates.LatticeConstant * [1; 1/sqrt(3)];
    b2 = 2*pi/cCoordinates.LatticeConstant * [0; 2/sqrt(3)];

    % special points in the BZ
    GammaPoint = [0;0];
    Kpoint_p = (2*b1-b2)/3;
    Kpoint_m = -(2*b1-b2)/3;
    Mpoint = b1/2;


    % creating vectors alon the Brillouine zone
    k_points_num = 100;
    k_section = zeros(2, 3*k_points_num);
    % Gamma -> M
    k_section(1, 1:k_points_num) = GammaPoint(1):(Mpoint(1)-GammaPoint(1))/(k_points_num-1):Mpoint(1);
    k_section(2, 1:k_points_num) = GammaPoint(2):(Mpoint(2)-GammaPoint(2))/(k_points_num-1):Mpoint(2);

    % M -> K
    k_section(1, k_points_num+1:2*k_points_num) = Mpoint(1):(Kpoint_p(1)-Mpoint(1))/(k_points_num-1):Kpoint_p(1);
    k_section(2, k_points_num+1:2*k_points_num) = Mpoint(2):(Kpoint_p(2)-Mpoint(2))/(k_points_num-1):Kpoint_p(2);

    % K -> Gamma
    k_section(1, 2*k_points_num+1:end) = Kpoint_p(1):(GammaPoint(1)-Kpoint_p(1))/(k_points_num-1):GammaPoint(1);
    k_section(2, 2*k_points_num+1:end) = zeros(1,k_points_num);
    
    % the length of the trajectory along the k_section
    k_length = zeros(1, 3*k_points_num);
    difference1 = norm([k_section(1, 2)-k_section(1, 1); k_section(2, 2)-k_section(2, 1)]);
    difference2 = norm([k_section(1, k_points_num+2)-k_section(1, k_points_num+1); k_section(2, k_points_num+2)-k_section(2, k_points_num+1)]);
    difference3 = norm([k_section(1, 2*k_points_num+2)-k_section(1, 2*k_points_num+1); k_section(2, 2*k_points_num+2)-k_section(2, 2*k_points_num+1)]);
    k_length(1:k_points_num) = (0:k_points_num-1)*difference1;
    k_length(k_points_num+1:2*k_points_num) = (0:k_points_num-1)*difference2 + k_points_num*difference1;
    k_length(2*k_points_num+1:3*k_points_num) = (0:k_points_num-1)*difference3 + k_points_num*(difference1+difference2);
    
    % delete duplicated points
    k_section(:,k_points_num) = [];
    k_section(:,2*k_points_num) = [];
    k_length(k_points_num) = [];
    k_length(2*k_points_num) = [];

    % calculate the spectra of monolayer TMDC (nearest and second neasrest neighbor models)
    [SpectralData1, SpectralData2] = SpectrumCalc();

    % plot the calculated data
    PlotFunction( )

    % exporting the calculated data
    save( fullfile(outputdir,[outfilename, '.mat']) );


%% CalcHamiltonian_k
%> @brief Calculates the spectrum on a triangle lattice from Fourier-transformed Hamiltonian
%> @param kvector a two-component column vector of the momentum vector in the BZ.
%> @return Returns with the calculated energy eigenvalue at a specific point determined by kvector.
    function ret = CalcHamiltonian_k( kvector )
        orbitals_num = 11;
        
        ret = zeros( orbitals_num );
        
        % obtaining derived hopping parameters t_2__i_i
        t_2__hoppings = param.scatter.Calc_t_2__i_i();
        
        % ******* EQ (4) in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a> *********
        for idx = 1:orbitals_num
    
            parameter_extension = ['__', num2str(idx), '_', num2str(idx)];
            epsilon = param.scatter.(['epsilon', num2str(idx)]);
            t_1 = param.scatter.(['t_1', parameter_extension]);
            t_2 = t_2__hoppings.(['t_2', parameter_extension]);
            ret( idx, idx) = ret( idx, idx) + epsilon + 2*t_1*cos(kvector'*delta1) + 2*t_2*(cos(kvector'*delta2) + cos(kvector'*delta3));
            
        end       
        
        
        % ******* EQ (5) in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a> *********
        % obtaining derived hopping parameters t_2__i_j
        t_2__hoppings = param.scatter.Calc_t_2__i_j();
        
        % obtaining derived hopping parameters t_3__i_j
        t_3__hoppings = param.scatter.Calc_t_3__i_j();
        
        % M_2 symmetry (+)
        orbital_indexes1 = [3, 6, 9];
        orbital_indexes2 = [5, 8, 11]; 
        for idx = 1:length( orbital_indexes1 )
            orbital_index1 = orbital_indexes1(idx);
            orbital_index2 = orbital_indexes2(idx);
            
            t_1 = param.scatter.(['t_1__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            t_2 = t_2__hoppings.(['t_2__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            t_3 = t_3__hoppings.(['t_3__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            ret(orbital_index1, orbital_index2) = ret(orbital_index1, orbital_index2) + 2*t_1*cos(kvector'*delta1) + ...
                            t_2*( exp(-1i*kvector'*delta2) + exp(-1i*kvector'*delta3) ) +...
                            t_3*( exp(1i*kvector'*delta2) + exp(1i*kvector'*delta3) );
                        
            ret(orbital_index2, orbital_index1) = ret(orbital_index2, orbital_index1) + 2*t_1*cos(kvector'*delta1) + ...
                            t_2*( exp(1i*kvector'*delta2) + exp(1i*kvector'*delta3) ) +...
                            t_3*( exp(-1i*kvector'*delta2) + exp(-1i*kvector'*delta3) );
              
        end
              
        
        
        % ******* EQ (6) in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a> *********
        % obtaining derived hopping parameters t_2__i_j
        t_2__hoppings = param.scatter.Calc_t_2__i_j();
        
        % obtaining derived hopping parameters t_3__i_j
        t_3__hoppings = param.scatter.Calc_t_3__i_j();
        
        % M_2 symmetry (-)
        orbital_indexes1 = [1, 3, 4, 6, 7, 9, 10];
        orbital_indexes2 = [2, 4, 5, 7, 8, 10, 11];  
        for idx = 1:length( orbital_indexes1 )
            orbital_index1 = orbital_indexes1(idx);
            orbital_index2 = orbital_indexes2(idx);
            
            t_1 = param.scatter.(['t_1__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            t_2 = t_2__hoppings.(['t_2__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            t_3 = t_3__hoppings.(['t_3__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            ret(orbital_index1, orbital_index2) = ret(orbital_index1, orbital_index2) - 2*1i*t_1*sin(kvector'*delta1) + ...
                            t_2*( exp(-1i*kvector'*delta2) - exp(-1i*kvector'*delta3) ) +...
                            t_3*( -exp(1i*kvector'*delta2) + exp(1i*kvector'*delta3) );
                        
            ret(orbital_index2, orbital_index1) = ret(orbital_index2, orbital_index1) + 2*1i*t_1*sin(kvector'*delta1) + ...
                            t_2*( exp(1i*kvector'*delta2) - exp(1i*kvector'*delta3) ) +...
                            t_3*( -exp(-1i*kvector'*delta2) + exp(-1i*kvector'*delta3) );
              
        end
        
        
        % ******* EQ (7) in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a> *********
        % obtaining derived hopping parameters t_4__i_j
        t_4__hoppings = param.scatter.Calc_t_4__i_j();
        
        
        % M_2 symmetry (+)
        orbital_indexes1 = [3, 5, 4, 10, 9, 11, 10];
        orbital_indexes2 = [1, 1, 2, 6, 7, 7, 8];  
        for idx = 1:length( orbital_indexes1 )
            orbital_index1 = orbital_indexes1(idx);
            orbital_index2 = orbital_indexes2(idx);
            
            t_4 = t_4__hoppings.(['t_4__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            ret(orbital_index1, orbital_index2) = ret(orbital_index1, orbital_index2) + t_4*( exp(1i*kvector'*delta4) - exp(1i*kvector'*delta6) );            
            ret(orbital_index2, orbital_index1) = ret(orbital_index2, orbital_index1) + t_4*( exp(-1i*kvector'*delta4) - exp(-1i*kvector'*delta6) );
              
        end
        
        
        % ******* EQ (8) in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a> *********        
        
        % M_2 symmetry (-)
        orbital_indexes1 = [4, 3, 5, 9, 11, 10, 9, 11];
        orbital_indexes2 = [1, 2, 2, 6, 6, 7, 8, 8];  
        for idx = 1:length( orbital_indexes1 )
            orbital_index1 = orbital_indexes1(idx);
            orbital_index2 = orbital_indexes2(idx);
            
            t_4 = t_4__hoppings.(['t_4__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            t_5 = param.scatter.(['t_5__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
            ret(orbital_index1, orbital_index2) = ret(orbital_index1, orbital_index2) + t_4*( exp(1i*kvector'*delta4) + exp(1i*kvector'*delta6) ) + ...
                                                    t_5*exp(1i*kvector'*delta5);             
            ret(orbital_index2, orbital_index1) = ret(orbital_index2, orbital_index1) + t_4*( exp(-1i*kvector'*delta4) + exp(-1i*kvector'*delta6) ) + ...
                                                    t_5*exp(-1i*kvector'*delta5);
              
        end   
        
        % ******* EQ (10) in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a> *********        
        ret(9,6) = ret(9,6) + param.scatter.('t_6__9_6')*(exp(1i*kvector'*delta7) + exp(1i*kvector'*delta8) + exp(1i*kvector'*delta9));
        ret(6,9) = ret(6,9) + param.scatter.('t_6__9_6')*(exp(-1i*kvector'*delta7) + exp(-1i*kvector'*delta8) + exp(-1i*kvector'*delta9));
        
        ret(11,6) = ret(11,6) + param.scatter.('t_6__11_6')*(exp(1i*kvector'*delta7) - 1/2*exp(1i*kvector'*delta8) - 1/2*exp(1i*kvector'*delta9));
        ret(6,11) = ret(6,11) + param.scatter.('t_6__11_6')*(exp(-1i*kvector'*delta7) - 1/2*exp(-1i*kvector'*delta8) - 1/2*exp(-1i*kvector'*delta9));
        
        ret(10,6) = ret(10,6) + param.scatter.('t_6__11_6')*sqrt(3)/2*(-exp(1i*kvector'*delta8) + exp(1i*kvector'*delta9));
        ret(6,10) = ret(6,10) + param.scatter.('t_6__11_6')*sqrt(3)/2*(-exp(-1i*kvector'*delta8) + exp(-1i*kvector'*delta9));
        
        ret(9,8) = ret(9,8) + param.scatter.('t_6__9_8')*(exp(1i*kvector'*delta7) - 1/2*exp(1i*kvector'*delta8) - 1/2*exp(1i*kvector'*delta9));
        ret(8,9) = ret(8,9) + param.scatter.('t_6__9_8')*(exp(-1i*kvector'*delta7) - 1/2*exp(-1i*kvector'*delta8) - 1/2*exp(-1i*kvector'*delta9));
        
        ret(9,7) = ret(9,7) + param.scatter.('t_6__9_8')*sqrt(3)/2*(-exp(1i*kvector'*delta8) + exp(1i*kvector'*delta9));
        ret(7,9) = ret(7,9) + param.scatter.('t_6__9_8')*sqrt(3)/2*(-exp(-1i*kvector'*delta8) + exp(-1i*kvector'*delta9));
        
        ret(10,7) = ret(10,7) + param.scatter.('t_6__11_8')*3/4*(exp(1i*kvector'*delta8) + exp(1i*kvector'*delta9));
        ret(7,10) = ret(7,10) + param.scatter.('t_6__11_8')*3/4*(exp(-1i*kvector'*delta8) + exp(-1i*kvector'*delta9));
        
        ret(11,7) = ret(11,7) + param.scatter.('t_6__11_8')*sqrt(3)/4*(exp(1i*kvector'*delta8) - exp(1i*kvector'*delta9));
        ret(7,11) = ret(7,11) + param.scatter.('t_6__11_8')*sqrt(3)/4*(exp(-1i*kvector'*delta8) - exp(-1i*kvector'*delta9));
        
        ret(10,8) = ret(10,8) + param.scatter.('t_6__11_8')*sqrt(3)/4*(exp(1i*kvector'*delta8) - exp(1i*kvector'*delta9));
        ret(8,10) = ret(8,10) + param.scatter.('t_6__11_8')*sqrt(3)/4*(exp(-1i*kvector'*delta8) - exp(-1i*kvector'*delta9));
        
        ret(11,8) = ret(11,8) + param.scatter.('t_6__11_8')*(exp(1i*kvector'*delta7) + 1/4*exp(1i*kvector'*delta8) + 1/4*exp(1i*kvector'*delta9));
        ret(8,11) = ret(8,11) + param.scatter.('t_6__11_8')*(exp(-1i*kvector'*delta7) + 1/4*exp(-1i*kvector'*delta8) + 1/4*exp(-1i*kvector'*delta9));
        
    end


%% SpectrumCalc
%> @brief Calculates the spectrum on a triangle lattice from Fourier-transformed Hamiltonian
%> @return [1] The calculated nearest neighbor spectrum.
%> @return [2] The calculated second nearest neighbor spectrum.
    function [SpectralData1, SpectralData2] = SpectrumCalc()
        
        % preallocating data arrays
        SpectralData1 = zeros( 11, size(k_section,2));
        SpectralData2 = zeros( 11, size(k_section,2));
        
        for idx = 1:size(k_section,2)
            
            % calculate the spectrum by the nearest neighbor model of EQuUs
            Hamiltonian = full(HLead.MomentumDependentHamiltonian( k_section(:,idx)'*a1, k_section(:,idx)'*a2));
            Eigenvalues = sort(real(eig(Hamiltonian)));
            SpectralData1(:, idx) = Eigenvalues;
            
            % calculate the spectrum by the second nearest neighbor model in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a>
            Hamiltonian = CalcHamiltonian_k( k_section(:,idx) );
            Eigenvalues = sort(real(eig(Hamiltonian)));
            SpectralData2(:, idx) = Eigenvalues;
        end

           
        
    end

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


 %% PlotFunction
%> @brief Creates the plot
    function PlotFunction( )
        
        % creating figure in units of pixels
        figure1 = figure( 'Units', 'Pixels');
        
        % font size on the figure will be 16 points
        fontsize = 16;
          
        
        
        axes_spectra = axes('Parent',figure1, ...
                'Visible', 'on',...
                'FontSize', fontsize,... 
                'Box', 'on',...
                'XTick', [], ...
                'Units', 'Pixels', ...
                'FontName','Times New Roman');
        hold on; 
        
        % plot the data (nearest neighbor model)
        for idx = 1:size(SpectralData1,1)
            EQuUs_plot = plot(k_length, SpectralData1(idx,:), 'Linewidth', 2, 'color', [0 0 0], 'Parent', axes_spectra); 
        end
        
        % plot the data (second nearest neighbor model)
        for idx = 1:size(SpectralData2,1)
            PRB_plot = plot(k_length, SpectralData2(idx,:), 'Linewidth', 1, 'color', [1 0 0], 'Parent', axes_spectra); 
        end
        
        
        % labeling special points in the BZ
        y_lim = get(axes_spectra, 'YLim');
        plot( k_length(k_points_num)*ones(1,2), y_lim, 'Linewidth', 1, 'color', [0 0 0], 'Parent', axes_spectra ) 
        plot( k_length(2*k_points_num-1)*ones(1,2), y_lim, 'Linewidth', 1, 'color', [0 0 0], 'Parent', axes_spectra ) 
        
        % latex interpreter is used only if display is detected
        if usejava('jvm') && feature('ShowFigureWindows')        
            text(0, -7, '$$\Gamma$$', 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_spectra, 'Interpreter', 'latex');
            text(k_length(k_points_num), -7, '$$M$$', 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_spectra, 'Interpreter', 'latex');
            text(k_length(2*k_points_num-1), -7, '$$K$$', 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_spectra, 'Interpreter', 'latex');
            text(k_length(3*k_points_num-10), -7, '$$\Gamma$$', 'FontSize', fontsize, 'FontName','Times New Roman',  'Parent', axes_spectra, 'Interpreter', 'latex');
        end
               
         
        
         % Create xlabel
         xlabel('k','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_spectra);

         % Create ylabel
         ylabel('Energy [eV]','FontSize', fontsize,'FontName','Times New Roman', 'Parent', axes_spectra);
         
         % set the legends
         fig_legend = legend( [EQuUs_plot, PRB_plot], {'EQuUs (nearest neighbor model)', 'PRB 92 205108 (second nearest neighbor model)'});%, 'FontSize', fontsize, 'FontName','Times New Roman')
         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(axes_spectra, 'Position', get(axes_spectra, 'OuterPosition') - get(axes_spectra, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1]);        
         Position_figure = get(axes_spectra, 'OuterPosition');
         set(figure1, 'Position', Position_figure);

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

end