test_TMDC_Bilayer_SOC.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_Bilayer_SOC.m
%> @brief Testfile to check the lattice construction of bilayer TMDC by class TMDC_Bilayer_SOC_Lead_Hamiltonians.
%> @Available
%> EQuUs v5.0 or later
%> @}
%> @brief Testfile to check the lattice construction of bilayer TMDC by class TMDC_Bilayer_SOC_Lead_Hamiltonians.
function test_TMDC_Bilayer_SOC( )


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

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

    % Parsing the input file and creating data structures
    [Opt, param] = parseInput( 'TMDC_Bilayer_Input_SOC.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;

    delta5 = (a1+2*a2)/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;

    
    % calculating the Hamiltonian at a special k point
    kvector = 0.3*Mpoint + 0.7*Kpoint_p;
    H_Fourier = CalcHamiltonian_k_Bilayer( kvector );
    
    % calculating the Hamiltonian at a special k point with EQuUs
    H_EQuUs = HLead.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);
    
    
    checkpoint = norm(  full(H_EQuUs) - full(H_Fourier) );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
    
    
    
%% CalcHamiltonian_k_Monolayer
%> @brief Calculates the spectrum on a TMDC monolayer lattice from Fourier-transformed Hamiltonian (PRB 92 205108)
%> @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_Monolayer( 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) );
            
            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.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);
                                                
            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   
        
        
        
        % ******  SOC: EQ (12) in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a> *********
        
        % adding spin degree of freedom
        ret  = kron( eye(2), ret);
        
        
        % defining orbital types according to Table II in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a>
        orbitals = {'d0_xz', 'd0_yz', 'p0_z', 'p0_x', 'p0_y', 'de_z2', 'de_xy', 'de_x2y2', 'pe_z', 'pe_x', 'pe_y'};
        orbital_number = length( orbitals );
        M = 1;
        site_number = orbital_number*M;
        SOC_H = zeros( site_number*2, site_number*2 );
        
        sigma_x = [0, 1; 1, 0]/2;
        sigma_y = [0, -1i; 1i, 0]/2;
        sigma_z = [1, 0; 0, -1]/2;
        
        for idx = 1:orbital_number
            orbital_type_row = orbitals{idx};
            row_indexes = (idx-1)*M+1:idx*M;
            
            for jdx = 1:orbital_number
                orbital_type_col = orbitals{jdx};
                col_indexes = (jdx-1)*M+1:jdx*M;
                
                % p and d orbitals are not coulped by SOC
                if ~strcmpi( orbital_type_row(1), orbital_type_col(1) )
                    continue
                end
                    
                
                if strcmpi( orbital_type_row(1), 'p' )
                % spin-orbit coupling of the p-type orbitals
                
                    if strcmpi( orbital_type_row, 'p0_z' ) && strcmpi( orbital_type_col, 'p0_x' )                        
                        SOC_H = SOC_H + kron( sigma_y, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 - ...
                                        kron( sigma_y, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;
                        
                    elseif strcmpi( orbital_type_row, 'p0_z' ) && strcmpi( orbital_type_col, 'p0_y' )                        
                        SOC_H = SOC_H + kron( sigma_x, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 - ...
                                        kron( sigma_x, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;
                                    
                    elseif strcmpi( orbital_type_row, 'p0_z' ) && strcmpi( orbital_type_col, 'pe_x' )                        
                        SOC_H = SOC_H + kron( sigma_y, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 + ...
                                        kron( sigma_y, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;
                                    
                    elseif strcmpi( orbital_type_row, 'p0_z' ) && strcmpi( orbital_type_col, 'pe_y' )                        
                        SOC_H = SOC_H + kron( sigma_x, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 + ...
                                        kron( sigma_x, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;   
                                    
                    elseif strcmpi( orbital_type_row, 'p0_x' ) && strcmpi( orbital_type_col, 'p0_y' )                        
                        SOC_H = SOC_H + kron( sigma_z, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 + ...
                                        kron( sigma_z, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;  
                                    
                    elseif strcmpi( orbital_type_row, 'p0_x' ) && strcmpi( orbital_type_col, 'pe_z' )                        
                        SOC_H = SOC_H + kron( sigma_y, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 + ...
                                        kron( sigma_y, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;  
                                    
                    elseif strcmpi( orbital_type_row, 'p0_x' ) && strcmpi( orbital_type_col, 'pe_y' )                        
                        SOC_H = SOC_H + kron( sigma_z, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 - ...
                                        kron( sigma_z, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;  
                        
                    elseif strcmpi( orbital_type_row, 'p0_y' ) && strcmpi( orbital_type_col, 'pe_z' )                        
                        SOC_H = SOC_H + kron( sigma_x, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 + ...
                                        kron( sigma_x, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;  
                                    
                    elseif strcmpi( orbital_type_row, 'p0_y' ) && strcmpi( orbital_type_col, 'pe_x' )                        
                        SOC_H = SOC_H + kron( sigma_z, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 - ...
                                        kron( sigma_z, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;
                                    
                    elseif strcmpi( orbital_type_row, 'pe_z' ) && strcmpi( orbital_type_col, 'pe_x' )                        
                        SOC_H = SOC_H + kron( sigma_y, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 - ...
                                        kron( sigma_y, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;
                                    
                    elseif strcmpi( orbital_type_row, 'pe_z' ) && strcmpi( orbital_type_col, 'pe_y' )                        
                        SOC_H = SOC_H + kron( sigma_x, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 - ...
                                        kron( sigma_x, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;
                                    
                    elseif strcmpi( orbital_type_row, 'pe_x' ) && strcmpi( orbital_type_col, 'pe_y' )                        
                        SOC_H = SOC_H + kron( sigma_z, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2 + ...
                                        kron( sigma_z, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_X, site_number, site_number) )/2;                                    
                                    
                    end
                
                
                elseif strcmpi( orbital_type_row(1), 'd' )
                % spin-orbit coupling of the d-type orbitals
                
                    if strcmpi( orbital_type_row, 'de_xy' ) && strcmpi( orbital_type_col, 'de_x2y2' )
                        SOC_H = SOC_H + kron( sigma_z, sparse( row_indexes, col_indexes, 2*1i*param.scatter.lambda_SO_M, site_number, site_number) );
                        
                    elseif strcmpi( orbital_type_row, 'de_xy' ) && strcmpi( orbital_type_col, 'd0_xz' )
                        SOC_H = SOC_H + kron( sigma_x, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_M, site_number, site_number) );
                        
                    elseif strcmpi( orbital_type_row, 'de_xy' ) && strcmpi( orbital_type_col, 'd0_yz' )
                        SOC_H = SOC_H + kron( sigma_y, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_M, site_number, site_number) );
                        
                    elseif strcmpi( orbital_type_row, 'de_x2y2' ) && strcmpi( orbital_type_col, 'd0_xz' )
                        SOC_H = SOC_H + kron( sigma_y, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_M, site_number, site_number) );
                        
                    elseif strcmpi( orbital_type_row, 'de_x2y2' ) && strcmpi( orbital_type_col, 'd0_yz' )
                        SOC_H = SOC_H + kron( sigma_x, sparse( row_indexes, col_indexes, 1i*param.scatter.lambda_SO_M, site_number, site_number) );
                        
                    elseif strcmpi( orbital_type_row, 'd0_xz' ) && strcmpi( orbital_type_col, 'd0_yz' )
                        SOC_H = SOC_H + kron( sigma_z, sparse( row_indexes, col_indexes, -1i*param.scatter.lambda_SO_M, site_number, site_number) );
                        
                    elseif strcmpi( orbital_type_row, 'd0_xz' ) && strcmpi( orbital_type_col, 'de_z2' )
                        SOC_H = SOC_H + kron( sigma_y, sparse( row_indexes, col_indexes, 1i*sqrt(3)*param.scatter.lambda_SO_M, site_number, site_number) );
                        
                    elseif strcmpi( orbital_type_row, 'd0_yz' ) && strcmpi( orbital_type_col, 'de_z2' )
                        SOC_H = SOC_H + kron( sigma_x, sparse( row_indexes, col_indexes, -1i*sqrt(3)*param.scatter.lambda_SO_M, site_number, site_number) );
                        
                    end
                    
                end
            
            end
        
        
        
        end
        
        SOC_H = SOC_H + SOC_H';
        
        ret = ret + SOC_H;
        
        
    end


%% CalcHamiltonian_k_Bilayer
%> @brief Calculates the spectrum on a TMDC bilayer lattice from Fourier-transformed Hamiltonian (PRB 92 205108)
%> @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_Bilayer( kvector )
        
        % upper layer
        % The Hamiltonian of the upper layer is tranformed by a mirror symmetry M3 in the xz plane (in this case mirroring the wave vector is equivalent to mirroring the lattice by y -> -y)
        % see the last paragraph of the left column on page 3 of PRB 92 205108
        upper_layer = CalcHamiltonian_k_Monolayer( [1,0;0,-1]*kvector );
        Umtx = diag( [1 -1 1 -1 1 1 -1 1 -1 1 -1] );
        Umtx_SOC = kron( [0, -1i; 1i, 0], Umtx ); % spin  operator should also be mirrored. The mirror operator to the xz plane is sigma_y        
        upper_layer = Umtx_SOC*upper_layer*Umtx_SOC';
        % lower layer
        lower_layer = CalcHamiltonian_k_Monolayer( kvector );
        %upper_layer = lower_layer;
        
        % ******  Bilyer EQ (12) in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a> *********
        
        ret  = [ upper_layer, zeros(size(lower_layer)); zeros(size(lower_layer)), lower_layer];
        
        
        % defining orbital types according to Table II in <a href="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.205108">PRB 92 205108</a>
        orbitals = {'d0_xz', 'd0_yz', 'p0_z', 'p0_x', 'p0_y', 'de_z2', 'de_xy', 'de_x2y2', 'pe_z', 'pe_x', 'pe_y'};
        
        orbital_number = length( orbitals );
        upper_site_num = length( orbitals );
        InterLayer_H = zeros( upper_site_num*2 );
        
        % relative shift between the upper and lower sheet
        stacking_vector = [delta5; (param.scatter.c-param.scatter.d_XX)]/param.scatter.a;
        
        
        
        for idx = 1:orbital_number
            % orbital type in the upper sheet
            orbital_type_upper = orbitals{idx};
            
            if strcmpi( orbital_type_upper(1), 'd' )
                % skipping the interlayer between the d orbitals 
                    continue
            end
            
            % site indexes of the upper layer
            upper_indexes = idx; 
            
            for jdx = 1:orbital_number
                % orbital type in the lower sheet
                orbital_type_lower = orbitals{jdx};                                    
                
                if strcmpi( orbital_type_lower(1), 'd' )
                    % skipping the interlayer between the d orbitals
                    continue
                end
                
                % site indexes of the lower layer
                lower_indexes = jdx + upper_site_num; 
                
                
                
                % ******************* interlayer hopping amplitudes within the unit cell *********************  
                
                % From lower sheet to upper sheet
                HoppingExtension = [ orbital_type_upper(end), orbital_type_lower(end) ];                
                r_vec = stacking_vector;    
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );                
                    InterLayer_H = InterLayer_H + sparse( upper_indexes, lower_indexes, InterlayerHoppings.(HoppingExtension)/2, upper_site_num*2, upper_site_num*2 );
                end
                
                
                % From upper sheet to lower sheet
                HoppingExtension = [ orbital_type_lower(end), orbital_type_upper(end) ];
                r_vec = -stacking_vector;     
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );                 
                    InterLayer_H = InterLayer_H + sparse( lower_indexes, upper_indexes, InterlayerHoppings.(HoppingExtension)/2, upper_site_num*2, upper_site_num*2 );
                end
                   
                
                                
                
                % ******************* interlayer hopping amplitudes between unit cells connected by a1 *********************
                
                % From upper sheet to lower sheet
                HoppingExtension = [ orbital_type_upper(end), orbital_type_lower(end) ];                
                r_vec = stacking_vector + [a1; 0]/param.scatter.a;   
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant ); 
                    InterLayer_H = InterLayer_H + sparse( lower_indexes, upper_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(1i*kvector'*a1), upper_site_num*2, upper_site_num*2 );
                end
                
                % From lower sheet to upper sheet
                HoppingExtension = [ orbital_type_lower(end), orbital_type_upper(end) ];  
                r_vec = -stacking_vector + [a1; 0]/param.scatter.a;   
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );
                    InterLayer_H = InterLayer_H + sparse( upper_indexes, lower_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(1i*kvector'*a1), upper_site_num*2, upper_site_num*2 );
                end
                
                % From upper sheet to lower sheet
                HoppingExtension = [ orbital_type_upper(end), orbital_type_lower(end) ];  
                r_vec = stacking_vector - [a1; 0]/param.scatter.a; 
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );                
                    InterLayer_H = InterLayer_H + sparse( lower_indexes, upper_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(-1i*kvector'*a1), upper_site_num*2, upper_site_num*2 );                                
                end
                
                % From lower sheet to upper sheet
                HoppingExtension = [ orbital_type_lower(end), orbital_type_upper(end) ];  
                r_vec = -stacking_vector - [a1; 0]/param.scatter.a;     
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );
                    InterLayer_H = InterLayer_H + sparse( upper_indexes, lower_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(-1i*kvector'*a1), upper_site_num*2, upper_site_num*2 );
                end
                
                
                                
                
                
                % ******************* interlayer hopping amplitudes between unit cells connected a2 *********************
                
                % From upper sheet to lower sheet
                HoppingExtension = [ orbital_type_upper(end), orbital_type_lower(end) ];                
                r_vec = stacking_vector + [a2; 0]/param.scatter.a;    
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant ); 
                    InterLayer_H = InterLayer_H + sparse( lower_indexes, upper_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(1i*kvector'*a2), upper_site_num*2, upper_site_num*2 );
                end
                
                % From lower sheet to upper sheet
                HoppingExtension = [ orbital_type_lower(end), orbital_type_upper(end) ];  
                r_vec = -stacking_vector + [a2; 0]/param.scatter.a;  
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );   
                    InterLayer_H = InterLayer_H + sparse( upper_indexes, lower_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(1i*kvector'*a2), upper_site_num*2, upper_site_num*2 );
                end
                
                % From upper sheet to lower sheet
                HoppingExtension = [ orbital_type_upper(end), orbital_type_lower(end) ];  
                r_vec = stacking_vector - [a2; 0]/param.scatter.a;
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );                
                    InterLayer_H = InterLayer_H + sparse( lower_indexes, upper_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(-1i*kvector'*a2), upper_site_num*2, upper_site_num*2 );
                end
                
                % From lower sheet to upper sheet
                HoppingExtension = [ orbital_type_lower(end), orbital_type_upper(end) ];  
                r_vec = -stacking_vector - [a2; 0]/param.scatter.a;  
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );
                    InterLayer_H = InterLayer_H + sparse( upper_indexes, lower_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(-1i*kvector'*a2), upper_site_num*2, upper_site_num*2 );
                end
                
                
                
                
                % ******************* interlayer hopping amplitudes between unit cells connected by a1+a2 *********************
                                
                % From upper sheet to lower sheet
                HoppingExtension = [ orbital_type_upper(end), orbital_type_lower(end) ];                
                r_vec = stacking_vector + [a1+a2; 0]/param.scatter.a;  
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );                
                    InterLayer_H = InterLayer_H + sparse( lower_indexes, upper_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(1i*kvector'*(a1+a2)), upper_site_num*2, upper_site_num*2 );
                end
                
                % From lower sheet to upper sheet
                HoppingExtension = [ orbital_type_lower(end), orbital_type_upper(end) ];  
                r_vec = -stacking_vector + [a1+a2; 0]/param.scatter.a;     
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );
                    InterLayer_H = InterLayer_H + sparse( upper_indexes, lower_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(1i*kvector'*(a1+a2)), upper_site_num*2, upper_site_num*2 );
                end
                
                % From upper sheet to lower sheet               
                HoppingExtension = [ orbital_type_upper(end), orbital_type_lower(end) ];  
                r_vec = stacking_vector - [a1+a2; 0]/param.scatter.a;   
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );                
                    InterLayer_H = InterLayer_H + sparse( lower_indexes, upper_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(-1i*kvector'*(a1+a2)), upper_site_num*2, upper_site_num*2 );
                end
                
                % From lower sheet to upper sheet
                HoppingExtension = [ orbital_type_lower(end), orbital_type_upper(end) ];  
                r_vec = -stacking_vector - [a1+a2; 0]/param.scatter.a; 
                if norm( r_vec*cCoordinates.LatticeConstant ) < 5
                    InterlayerHoppings = param.scatter.Calc_Hopping( r_vec*cCoordinates.LatticeConstant );
                    InterLayer_H = InterLayer_H + sparse( upper_indexes, lower_indexes, InterlayerHoppings.(HoppingExtension)/2*exp(-1i*kvector'*(a1+a2)), upper_site_num*2, upper_site_num*2 );
                end
                
                
                    
            end
        
        
        
        end
        
        % includng spin degree of freedoms
        InterLayer_H_SOC = sparse( [], [], [], upper_site_num*4, upper_site_num*4 );
        InterLayer_H_SOC( 1:upper_site_num, upper_site_num*2+1:upper_site_num*3 ) = InterLayer_H( 1:upper_site_num, upper_site_num+1:upper_site_num*2 );        
        InterLayer_H_SOC( upper_site_num+1:upper_site_num*2, upper_site_num*3+1:upper_site_num*4 ) = InterLayer_H( 1:upper_site_num, upper_site_num+1:upper_site_num*2 );
        
        InterLayer_H_SOC( upper_site_num*2+1:upper_site_num*3, 1:upper_site_num ) = InterLayer_H( upper_site_num+1:upper_site_num*2, 1:upper_site_num );
        InterLayer_H_SOC( upper_site_num*3+1:upper_site_num*4, upper_site_num+1:upper_site_num*2  ) = InterLayer_H( upper_site_num+1:upper_site_num*2, 1:upper_site_num );
        
       
        % adding interlayer coupling to the Hmailtonian
        ret = ret + InterLayer_H_SOC;
        
    end


end