test_TMDC_Monolayer.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_TMDC_Monolayer.m
%> @brief Testfile to check the lattice construction of monolayer TMDC by classes TMDC_Monolayer_Lead_Hamiltonians and TMDC_Monolayer_SOC_Lead_Hamiltonians.
%> @Available
%> EQuUs v5.0 or later
%> @}
%> @brief Testfile to check the lattice construction of monolayer TMDC by classes TMDC_Monolayer_Lead_Hamiltonians and TMDC_Monolayer_SOC_Lead_Hamiltonians.
function test_TMDC_Monolayer()


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

%% test of the class TMDC_Monolayer_Lead_Hamiltonians
    disp( '***** test of the class TMDC_Monolayer_Lead_Hamiltonians *****' );

    % create input structures
    [Opt, param] = parseInput( 'TMDC_Input.xml');

    % setting the width of the lead 1
    param.Leads{1}.M = 1;

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

    % getting lattice vectors of the triangle lattice
    cCoordinates = HLead_1.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;

    % Calculate the energy eigenvalue at a specific K-point
    kvector = [0.5; 0.3];
    %using EQuUs
    H_EQuUs = HLead_1.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);

    % using the Fourier-transformed Hamiltonian in PRB 92 205108
    H_Fourier = CalcHamiltonian_k( kvector );
    
    checkpoint = norm(  full(H_EQuUs) - full(H_Fourier) );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
%% test of the energy eigenvalues for lead of width M>1
    disp( '***** test of the energy eigenvalues for lead of width M>1 *****' );
    
    % create input structures
    [Opt, param] = parseInput( 'TMDC_Input.xml');

    % setting the width of the first and second leads
    param.Leads{1}.M = 1;
    param.Leads{2}.M = 3;
    
    % Calculate the energy eigenvalue at a specific K-point
    kvector = [0.5; 0.3];
    

    % creating class CreateLeadHamiltonians
    HLead_2 = CreateLeadHamiltonians(Opt, param, 'Hanyadik_Lead',2);
    
    % creating Hamiltonians
    HLead_2.CreateHamiltonians();
    
    % Creating the Hamiltonian of the M width lead.
    Hamiltonian_M3 = HLead_2.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2*param.Leads{2}.M);
    
    % Creating the Hamiltonian of the M=1 width lead.
    Hamiltonian_M1 = HLead_1.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);
    
    % calculate the eigenvector of Hamiltonai M1
    eigenvalues_M1 = real(eig(full(Hamiltonian_M1)));
    
    for jdx = 1:11
        eigval = eigenvalues_M1(jdx);
    
        difference = eig(full(Hamiltonian_M3)) - ones(size(Hamiltonian_M3,1),1)*eigval;
        difference = sort( abs(difference), 'ascend' );
    
        checkpoint = abs( difference(1) );
        if checkpoint > 1e-10
            warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
        end   
    end
    
    
    
%% test of the class TMDC_Monolayer_SOC_Lead_Hamiltonians
    disp( '***** test of the class TMDC_Monolayer_SOC_Lead_Hamiltonians *****' );  
    
    % create input structures
    [Opt, param] = parseInput( 'TMDC_Input_SOC.xml');

    % setting the width of the lead 1
    param.Leads{1}.M = 1;

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

    % getting lattice vectors of the triangle lattice
    cCoordinates = HLead_1.Read('coordinates');
    a1 = cCoordinates.a*cCoordinates.LatticeConstant;
    a2 = cCoordinates.b*cCoordinates.LatticeConstant;
    
    % 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]; 
    Mpoint = b1/2;

    % Calculate the Hamiltonian at a specific K-point
    kvector = 0.36*(Mpoint - GammaPoint);
    %using EQuUs
    H_EQuUs = HLead_1.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);

    % The spin splitting along the Gamma - M direction is zero.
    EigenValues = sort(real(eig(full(H_EQuUs))));
    DiffEigenValues = diff(EigenValues);
    DiffEigenValues = DiffEigenValues(1:2:end);
    
    checkpoint = norm(  DiffEigenValues );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end 
    
    
    
%% test of the energy eigenvalues for lead of width M>1
    disp( '***** test of the energy eigenvalues for lead of width M>1 *****' );
    
    % create input structures
    [Opt, param] = parseInput( 'TMDC_Input_SOC.xml');

    % setting the width of the first and second leads
    param.Leads{1}.M = 1;
    param.Leads{2}.M = 3;
    
    % Calculate the energy eigenvalue at a specific K-point
    kvector = [0.5; 0.3];
    

    % creating class CreateLeadHamiltonians
    HLead_2 = CreateLeadHamiltonians(Opt, param, 'Hanyadik_Lead',2);
    
    % creating Hamiltonians
    HLead_2.CreateHamiltonians();
    
    % Creating the Hamiltonian of the M width lead.
    Hamiltonian_M3 = HLead_2.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2*param.Leads{2}.M);
    
    % Creating the Hamiltonian of the M=1 width lead.
    Hamiltonian_M1 = HLead_1.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);
    
    % calculate the eigenvector of Hamiltonai M1
    eigenvalues_M1 = real(eig(full(Hamiltonian_M1)));
    
    for jdx = 1:22
        eigval = eigenvalues_M1(jdx);
    
        difference = eig(full(Hamiltonian_M3)) - ones(size(Hamiltonian_M3,1),1)*eigval;
        difference = sort( abs(difference), 'ascend' );
    
        checkpoint = abs( difference(1) );
        if checkpoint > 1e-10
            warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
        end   
    end    
    
    
    
    
    



%% 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.Leads{1}.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.Leads{1}.(['epsilon', num2str(idx)]);
            t_1 = param.Leads{1}.(['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.Leads{1}.Calc_t_2__i_j();
        
        % obtaining derived hopping parameters t_3__i_j
        t_3__hoppings = param.Leads{1}.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.Leads{1}.(['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.Leads{1}.Calc_t_2__i_j();
        
        % obtaining derived hopping parameters t_3__i_j
        t_3__hoppings = param.Leads{1}.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.Leads{1}.(['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.Leads{1}.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*( 1 - exp(1i*kvector'*delta1) );            
            ret(orbital_index2, orbital_index1) = ret(orbital_index2, orbital_index1) + t_4*( 1 - exp(-1i*kvector'*delta1) );
              
        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.Leads{1}.(['t_5__', num2str(orbital_index1), '_', num2str(orbital_index2)]);
                                                
            ret(orbital_index1, orbital_index2) = ret(orbital_index1, orbital_index2) + t_4*( 1 + exp(1i*kvector'*delta1) ) + ...
                                                    t_5*exp(1i*kvector'*delta2);            
            ret(orbital_index2, orbital_index1) = ret(orbital_index2, orbital_index1) + t_4*( 1 + exp(-1i*kvector'*delta1) ) + ...
                                                    t_5*exp(-1i*kvector'*delta2);
              
        end           
       
        
    end



end