Hex_Lead_Hamiltonians.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - Hex_Lead_Hamiltonians
%    Copyright (C) 2009-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 lattices Lattices
%> @{
%> @file Hex_Lead_Hamiltonians.m
%> @brief A class to create the Hamiltonian of one unit cell in a translational invariant lead made of hexagonal lattice structure, including graphene and silicene.
%> @} 
%> @brief A class to create the Hamiltonian of one unit cell in a translational invariant lead made of hexagonal lattice structure, including graphene and silicene.
%%
classdef Hex_Lead_Hamiltonians
    

methods ( Access = public )
    
 

%% Graphene_Bilayer_Lead_Hamiltonians
%> @brief Creates Hamiltonians H_0 and H_1 of bilayer graphene ribbon with zigzag/armchair edge as well as the structure conatining the coordinates of the atomic sites.
%> @image html bilayer_type1.jpg
%> @image latex bilayer_type1.jpg
%> @param lead_param An instance of structure #Lattice_Bilayer_Graphene (or its subclass) containing the physical parameters.
%> @param M Number of sites in the cross section of the lead.
%> @param End_Type The orientation of the lattice. Set 'A' for lattice with armchair orientation (meaning zigzag edges) or 'Z' for zizgaz orientation (meaning arm-chair edges)
%> @return [1] The Hamiltonian of one slab in the ribbon. 
%> @return [2] The coupling between the slabs. 
%> @return [3] The transverse coupling between the slabs for transverse calculations. 
%> @return [4] A structure Coordinates containing the coordinates of the sites. 
function [H0, H1, H1_transverse, coordinates] = Graphene_Bilayer_Lead_Hamiltonians(obj, lead_param, M, End_Type )
    
    % check the structure containing the physical parameters
    supClasses = superclasses(lead_param);
    if sum( strcmp( supClasses, 'Lattice_Graphene_Bilayer' ) ) == 0
        error(['EQuUs:Lattices:', class(obj), ':Graphene_Bilayer_Lead_Hamiltonians'], 'Invalid type of the input parameter');   
    end
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Graphene_Bilayer_Lead_Hamiltonians'], 'The input parameter M is empty')
    end    
    
    [H0_lower,H1_lower,H1_transverse_lower,coordinates_lower] = obj.Graphene_Lead_Hamiltonians(lead_param, M, End_Type);
    
    coordinates_lower.z = zeros(size(coordinates_lower.x));
    
    H0_upper = H0_lower;
    H1_upper = H1_lower;
    H1_transverse_upper = H1_transverse_lower;
    
    if isempty( lead_param.vargamma1 )
        lead_param.vargamma1 = 0.381;
    end
    
    if isempty( lead_param.vargamma3 )
        lead_param.vargamma3 = 0;
    end
    
    % zigzag edged bilayer ribbon: warping checked
    if strcmp(End_Type, 'A')
        Bernal_stack = [0, 1;-1,0]*coordinates_lower.a/norm(coordinates_lower.a); 
        
        H0_gamma1 = sparse(2:2:M, 1:2:M-1, -lead_param.vargamma1, size(H0_lower,1), size(H0_upper,2)) + sparse( M+3:2:size(H0_lower,1), M+2:2:size(H0_upper,2)-1, -lead_param.vargamma1, size(H0_lower,1), size(H0_upper,2)); 
        H0_gamma3 = sparse(1:2:M-1, M+1:2:size(H0_upper,2)-1, -lead_param.vargamma3, size(H0_lower,1), size(H0_upper,2)) + ... 
            sparse(M+2:2:size(H0_lower,1), M+1:2:size(H0_upper,2)-1, -lead_param.vargamma3, size(H0_lower,1), size(H0_upper,2)) + ...
            sparse(M+2:2:size(H0_lower,1), 2:2:M, -lead_param.vargamma3, size(H0_lower,1), size(H0_upper,2)) + ...
            sparse(3:2:M-1, 2:2:M-1, -lead_param.vargamma3, size(H0_lower,1), size(H0_upper,2));
        H0 = [H0_lower, H0_gamma1 + H0_gamma3;...
            (H0_gamma1+H0_gamma3)' , H0_upper];
    
        H1_gamma3_lower_upper = sparse(M+2:2:size(H1_lower,1), 2:2:M, -lead_param.vargamma3, size(H1_lower,1), size(H1_upper,2));
        H1_gamma3_upper_lower = sparse(M+1:2:size(H1_upper,1)-1, 1:2:M-1, -lead_param.vargamma3, size(H1_upper,1), size(H1_lower,2));
        H1 = [H1_lower, H1_gamma3_lower_upper; ...
                H1_gamma3_upper_lower, H1_upper];
            
          
        if ~isempty(H1_transverse_lower)
            %> tranverse coupling: warping checked
            H1_transverse_lower_upper = sparse(M+1, 2*M, -lead_param.vargamma1, size(H1_transverse_lower,1), size(H1_transverse_upper,2));
            H1_transverse_upper_lower = sparse(2*M, M+1, -lead_param.vargamma1, size(H1_transverse_upper,1), size(H1_transverse_lower,2));
            H1_transverse_gamma3_lower_upper = sparse(1, M, -lead_param.vargamma3, size(H1_transverse_lower,1), size(H1_transverse_upper,2));
            H1_transverse_gamma3_upper_lower = sparse([], [], [], size(H1_transverse_upper,1), size(H1_transverse_lower,2));
            H1_transverse = [ H1_transverse_lower, H1_transverse_lower_upper + H1_transverse_gamma3_lower_upper;
                        H1_transverse_upper_lower + H1_transverse_gamma3_upper_lower, H1_transverse_upper];
        else
            H1_transverse = [];
        end
            
            
    % armchair edged bilayer ribbon: warping checked, but not for H1_transverse
    elseif strcmp(End_Type, 'Z')
        
        Bernal_stack = [-0.5; sqrt(3)/2]; 
        
        H0_gamma1 = sparse(1:M, M+1:2*M, -lead_param.vargamma1, size(H0_lower,1), size(H0_upper,2)) + sparse( 2*M+2:3*M, 3*M+1:4*M-1, -lead_param.vargamma1, size(H0_lower,1), size(H0_upper,2)); 
        H0_gamma3 = sparse(M+2:2*M, 1:M-1, -lead_param.vargamma3, size(H0_lower,1), size(H0_upper,2)) + ...
            sparse(M+1:2*M, 2*M+1:3*M, -lead_param.vargamma3, size(H0_lower,1), size(H0_upper,2)) + sparse(M+2:2*M, 2*M+1:3*M-1, -lead_param.vargamma3, size(H0_lower,1), size(H0_upper,2)) + ...
            sparse(3*M+1:4*M, 2*M+1:3*M, -lead_param.vargamma3, size(H0_lower,1), size(H0_upper,2));
        H0 = [H0_lower, H0_gamma1 + H0_gamma3;...
            (H0_gamma1+H0_gamma3)' , H0_upper];
    
        H1_gamma3_lower_upper = sparse(3*M+1:4*M, 1:M, -lead_param.vargamma3, size(H1_lower,1), size(H1_upper,2)) + sparse(3*M+2:4*M, 1:M-1, -lead_param.vargamma3, size(H1_lower,1), size(H1_upper,2));
        H1_gamma3_upper_lower = sparse([], [], -lead_param.vargamma3, size(H1_upper,1), size(H1_lower,2));
        H1 = [H1_lower, H1_gamma3_lower_upper; ...
                H1_gamma3_upper_lower, H1_upper];
            
            
            
        if ~isempty(H1_transverse_lower)
            % warping needs to be included
            H1_transverse_lower_upper = sparse(2*M+1, 4*M, -lead_param.vargamma1, size(H1_transverse_lower,1), size(H1_transverse_upper,2));
            H1_transverse_upper_lower = sparse(4*M, 2*M+1, -lead_param.vargamma1, size(H1_transverse_upper,1), size(H1_transverse_lower,2));
            H1_transverse = [ H1_transverse_lower, H1_transverse_lower_upper;
                        H1_transverse_upper_lower, H1_transverse_upper];
        else
            H1_transverse = [];
        end
        
        
    else
        save('Hex_Lead_Hamiltonians_Graphene_Bilayer_Lead_Hamiltonians.mat');
        error('Undefined lead end type')
    end
    
                   
    coordinates_upper = structures( 'coordinates' );
    fields = fieldnames( coordinates_lower );
    for idx = 1:length( fields )
        fieldname = fields{idx};
        if strcmp( fieldname, 'x' )
            coordinates_upper.x = coordinates_lower.x + Bernal_stack(1);            
        elseif strcmp( fieldname, 'y' )            
            coordinates_upper.y = coordinates_lower.y + Bernal_stack(2);
        elseif strcmp( fieldname, 'z' )            
            coordinates_upper.z = coordinates_lower.z + 3.3/1.42;    %The distance between the layers in units of rCC is 3.3 \AA ?
        else
            coordinates_upper.(fieldname) = coordinates_lower.(fieldname);
        end
    end
    

    coordinates = structures( 'coordinates' );  
    for idx = 1:length( fields )
        fieldname = fields{idx};
        if strcmp( fieldname, 'a' )
            coordinates.a  = coordinates_lower.(fieldname);
        elseif strcmp( fieldname, 'b' )
            coordinates.b  = coordinates_lower.(fieldname);
        elseif strcmp( fieldname, 'LatticeConstant' )
            coordinates.LatticeConstant  = coordinates_lower.LatticeConstant;
        else
            coordinates.(fieldname) = [coordinates_lower.(fieldname); coordinates_upper.(fieldname)];
        end
    end
    
end



%% Graphene_Bilayer_Lead_Hamiltonians2
%> @brief Creates Hamiltonians H_0 and H_1 of bilayer graphene ribbon type 2 with zigzag/armchair edge as well as the structure conatining the coordinates of the atomic sites.
%> @image html bilayer_type2.jpg
%> @image latex bilayer_type2.jpg
%> @param lead_param An instance of structure #Lattice_Bilayer_Graphene (or its subclass) containing the physical parameters.
%> @param M Number of sites in the cross section of the lead.
%> @param End_Type The orientation of the lattice. Set 'A' for lattice with armchair orientation (meaning zigzag edges) or 'Z' for zizgaz orientation (meaning arm-chair edges)
%> @return [1] The Hamiltonian of one slab in the ribbon. 
%> @return [2] The coupling between the slabs. 
%> @return [3] The transverse coupling between the slabs for transverse calculations. 
%> @return [4] A structure Coordinates containing the coordinates of the sites. 
function [H0, H1, H1_transverse, coordinates] = Graphene_Bilayer_Lead_Hamiltonians2(obj, lead_param, M, End_Type, varargin )
    
    % check the structure containing the physical parameters
    supClasses = superclasses(lead_param);
    if sum( strcmp( supClasses, 'Lattice_Graphene_Bilayer' ) ) == 0
        error(['EQuUs:Lattices:', class(obj), ':Graphene_Bilayer_Lead_Hamiltonians2'], 'Invalid type of the input parameter');   
    end
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Graphene_Bilayer_Lead_Hamiltonians2'], 'The input parameter M is empty')
    end
        
    
    % zigzag edged bilayer ribbon
    if strcmp(End_Type, 'A')
    
        Mtmp = M+4;
    
        [H0, H1, H1_transverse, coordinates] = obj.Graphene_Bilayer_Lead_Hamiltonians(lead_param, Mtmp, End_Type, varargin{:});
    
        sites_to_be_removed = [1, 2, Mtmp+1, Mtmp+2, Mtmp-1, Mtmp, 2*Mtmp, ... lower sheet
                                        [1, 2, Mtmp+1, Mtmp-1, Mtmp, 2*Mtmp-1, 2*Mtmp] + 2*Mtmp]; %upper sheet
                                    
        logical_indexes = true(size(H0, 1),1);
        logical_indexes(sites_to_be_removed) = false;
    
        H0 = H0(logical_indexes, logical_indexes);
        H1 = H1(logical_indexes, logical_indexes);
    
        if ~isempty(H1_transverse)
            H1_transverse = H1_transverse(logical_indexes, logical_indexes);
        end
    
    
        fields = fieldnames( coordinates );
        for idx = 1:length( fields )
            fieldname = fields{idx};
            if strcmpi( fieldname, 'a') || strcmpi( fieldname, 'b') || isempty(coordinates.(fieldname)) || strcmpi( fieldname, 'LatticeConstant')
                continue
            end
        
            coordinates.(fieldname) = coordinates.(fieldname)(logical_indexes);
        
        end
    
    
    % armchair edged bilayer ribbon
    elseif strcmp(End_Type, 'Z')
        [H0, H1, H1_transverse, coordinates] = obj.Graphene_Bilayer_Lead_Hamiltonians(lead_param, M, End_Type, varargin{:});        
    end
    
    
end



%% Graphene_Bilayer_Lead_Hamiltonians3
%> @brief Creates Hamiltonians H_0 and H_1 of bilayer graphene ribbon with zigzag/armchair edge as well as the structure conatining the coordinates of the atomic sites.
%> @param lead_param An instance of structure #Lattice_Bilayer_Graphene (or its subclass) containing the physical parameters.
%> @param M Number of sites in the cross section of the lead.
%> @param End_Type The orientation of the lattice. Set 'A' for lattice with armchair orientation (meaning zigzag edges) or 'Z' for zizgaz orientation (meaning arm-chair edges)
%> @param varargin Optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'q' The transverse momentum. Set to empty (default) for computations without transverse momentums.
%> @return [1] The Hamiltonian of one slab in the ribbon. 
%> @return [2] The coupling between the slabs. 
%> @return [3] The transverse coupling between the slabs for transverse calculations. 
%> @return [4] A structure Coordinates containing the coordinates of the sites. 
function [H0, H1, H1_transverse, coordinates] = Graphene_Bilayer_Lead_Hamiltonians3(obj, lead_param, M, End_Type )
    
    % check the structure containing the physical parameters
    supClasses = superclasses(lead_param);
    if sum( strcmp( supClasses, 'Lattice_Graphene_Bilayer' ) ) == 0
        error(['EQuUs:Lattices:', class(obj), ':Graphene_Bilayer_Lead_Hamiltonians3'], 'Invalid type of the input parameter');   
    end
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Graphene_Bilayer_Lead_Hamiltonians3'], 'The input parameter M is empty')
    end
        
    Mtmp = M+2;
    
    [H0, H1, H1_transverse, coordinates] = obj.Graphene_Bilayer_Lead_Hamiltonians(lead_param, Mtmp, End_Type );
    
    sites_to_be_removed = [1, Mtmp+1, ... lower sheet
                            [Mtmp, 2*Mtmp] + 2*Mtmp]; %upper sheet
                                    
    logical_indexes = true(size(H0, 1),1);
    logical_indexes(sites_to_be_removed) = false;
    
    H0 = H0(logical_indexes, logical_indexes);
    H1 = H1(logical_indexes, logical_indexes);
    
    if ~isempty(H1_transverse)
        H1_transverse = H1_transverse(logical_indexes, logical_indexes);
    end
    
    
    fields = fieldnames( coordinates );
    for idx = 1:length( fields )
        fieldname = fields{idx};
        if strcmpi( fieldname, 'a') || strcmpi( fieldname, 'b') || isempty(coordinates.(fieldname))
            continue
        end
        
        coordinates.(fieldname) = coordinates.(fieldname)(logical_indexes);
        
    end
    
    
    
end

%% Silicene_Lead_Hamiltonians
%> @brief Creates Hamiltonians H_0 and H_1 of silicene ribbon with zigzag edge as well as the structure conatining the coordinates of the atomic sites.
%> @param lead_param An instance of structure #Lattice_Silicene (or its subclass) containing the physical parameters.
%> @param M Number of sites in the cross section of the lead.
%> @param End_Type The orientation of the lattice. Set 'A' for lattice with armchair orientation (meaning zigzag edges) or 'Z' for zizgaz orientation (meaning arm-chair edges)
%> @return [1] The Hamiltonian of one slab in the ribbon. 
%> @return [2] The coupling between the slabs. 
%> @return [3] The transverse coupling between the slabs for transverse calculations. 
%> @return [4] A structure Coordinates containing the coordinates of the sites. 
function [H0, H1, H1_transverse, coordinates] = Silicene_Lead_Hamiltonians(obj, lead_param, M, End_Type )
    % Hamiltonian created according to arXiv:1305.3684v1
    
    % check the structure containing the physical parameters
    supClasses = superclasses(lead_param);
    if sum( strcmp( supClasses, 'Lattice_Silicene' ) ) == 0
        error(['EQuUs:Lattices:', class(obj), ':Silicene_Lead_Hamiltonians'], 'Invalid type of the input parameter');   
    end
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Silicene_Lead_Hamiltonians'], 'The input parameter M is empty')
    end    
    
    if strcmp( End_Type, 'Z' )
        error(['EQuUs:', class(obj), ':Silicene_Lead_Hamiltonians'], 'Silicene Hamiltonian developed only for Armchair ended leads')
    end
    

    if isprop( lead_param, 'SOintrinsic' )
        SOintrinsic = -lead_param.SOintrinsic/(3*sqrt(3));
    else
        SOintrinsic = -0.3*vargamma/(3*sqrt(3));
    end
    if isprop( lead_param, 'SO_rashba_intrinsic' )
        SO_rashba_intrinsic = lead_param.SO_rashba_intrinsic*2/3;
    else
        SO_rashba_intrinsic = 0.07*vargamma;
    end
    
    
    sgn = 1;

    [H0,H1,H1_transverse,coordinates] = Graphene_Lead_Hamiltonians(obj, lead_param, M, End_Type);
    coordinates.spinup = [ true( size(coordinates.x) ); false( size(coordinates.x) ) ];
    coordinates.x = [ coordinates.x; coordinates.x];
    coordinates.y = [ coordinates.y; coordinates.y];
    
    
    % intrinsic SO coupling
    HSOintrinsic_H0 = sparse( 1:2:M-1, M+2:2:2*M, -1i*SOintrinsic, size(H0,1), size(H0,2)) + sparse( 2:2:M, M+1:2:2*M-1, 1i*SOintrinsic, size(H0,1), size(H0,2)) + ...
                   sparse( 3:2:M, M+2:2:2*M-1, -1i*SOintrinsic, size(H0,1), size(H0,2)) + sparse( 2:2:M-1, M+3:2:2*M, 1i*SOintrinsic, size(H0,1), size(H0,2));
    HSOintrinsic_H0 = HSOintrinsic_H0 + HSOintrinsic_H0';
    
    HSOintrinsic_H1 = sparse( 1:2:M, 1:2:M, 1i*SOintrinsic, size(H1,1), size(H1,2)) + sparse( 2:2:M, 2:2:M, -1i*SOintrinsic, size(H1,1), size(H1,2)) + ...
                   sparse( M+1:2:2*M, M+1:2:2*M, -1i*SOintrinsic, size(H1,1), size(H1,2)) + sparse( M+2:2:2*M, M+2:2:2*M, 1i*SOintrinsic, size(H1,1), size(H1,2)) + ...
                   sparse( M+1:2:2*M, 2:2:M, 1i*SOintrinsic, size(H1,1), size(H1,2)) + sparse( M+2:2:2*M-1, 3:2:M, -1i*SOintrinsic, size(H1,1), size(H1,2)) + ...
                   sparse( M+3:2:2*M, 2:2:M-1, 1i*SOintrinsic, size(H1,1), size(H1,2)) + sparse( M+2:2:2*M, 1:2:M, -1i*SOintrinsic, size(H1,1), size(H1,2));
               
              
    % intrinsic Rashba SO coupling
    % next-nearest vectors
    d1 = [coordinates.x(M+2) - coordinates.x(1); coordinates.y(M+2) - coordinates.y(1)]; d1 = d1/norm(d1);
    d2 = [coordinates.x(M+1) - coordinates.x(2); coordinates.y(M+1) - coordinates.y(2)]; d2 = d2/norm(d2);
    d3 = coordinates.a; d3 = d3/norm(d3);
    
    muvec = 0:M-1;
    muvec_firstrow = (-1).^(muvec) * sgn;
    muvec_secondrow = -muvec_firstrow;
    % matrix of spin updown
    HSO_rashba_intrinsic_H0 = ...
        sparse( M+2:2*M, 1:M-1, -1i*SO_rashba_intrinsic*(d1(2)+1i*d1(1))*muvec_secondrow(2:M), size(H0,1), size(H0,2)) + ...
        sparse( 1:M-1, M+2:2*M, -1i*SO_rashba_intrinsic*(-d1(2)+1i*(-d1(1)))*(muvec_firstrow(1:M-1)), size(H0,1), size(H0,2)) + ...
        sparse( M+1:2*M-1, 2:M, -1i*SO_rashba_intrinsic*(d2(2)+1i*d2(1))*muvec_secondrow(1:M-1), size(H0,1), size(H0,2)) + ...
        sparse( 2:M, M+1:2*M-1, -1i*SO_rashba_intrinsic*(-d2(2)+1i*(-d2(1)))*(muvec_firstrow(2:M)), size(H0,1), size(H0,2));
    HSO_rashba_intrinsic_H1 = ...
        sparse( 1:M, 1:M,         -1i*SO_rashba_intrinsic*(d3(2)+1i*d3(1))*muvec_firstrow, size(H1,1), size(H1,2)) + ...
        sparse( M+1:2*M, M+1:2*M, -1i*SO_rashba_intrinsic*(d3(2)+1i*d3(1))*muvec_secondrow, size(H1,1), size(H1,2)) + ...
        sparse( 2:M, M+1:2*M-1,         -1i*SO_rashba_intrinsic*(d1(2)+1i*d1(1))*muvec_firstrow(2:M), size(H1,1), size(H1,2)) + ...
        sparse( 1:M-1, M+2:2*M,         -1i*SO_rashba_intrinsic*(d2(2)+1i*d2(1))*muvec_firstrow(1:M-1), size(H1,1), size(H1,2));
    
    % matrix of spin downup    
    HSO_rashba_intrinsic_H0b = ...
        sparse( M+2:2*M, 1:M-1, -1i*SO_rashba_intrinsic*(d1(2)-1i*d1(1))*muvec_secondrow(2:M), size(H0,1), size(H0,2)) + ...
        sparse( 1:M-1, M+2:2*M, -1i*SO_rashba_intrinsic*(-d1(2)-1i*(-d1(1)))*muvec_firstrow(1:M-1), size(H0,1), size(H0,2)) + ...
        sparse( M+1:2*M-1, 2:M, -1i*SO_rashba_intrinsic*(d2(2)-1i*d2(1))*muvec_secondrow(1:M-1), size(H1,1), size(H0,2)) + ...
        sparse( 2:M, M+1:2*M-1, -1i*SO_rashba_intrinsic*(-d2(2)-1i*(-d2(1)))*muvec_firstrow(2:M), size(H1,1), size(H0,2));
    HSO_rashba_intrinsic_H1b = ...
        sparse( 1:M, 1:M,         -1i*SO_rashba_intrinsic*(d3(2)-1i*d3(1))*muvec_firstrow, size(H0,1), size(H1,2)) + ...
        sparse( M+1:2*M, M+1:2*M, -1i*SO_rashba_intrinsic*(d3(2)-1i*d3(1))*muvec_secondrow, size(H0,1), size(H1,2)); 
        sparse( 2:M, M+1:2*M-1,         -1i*SO_rashba_intrinsic*(d1(2)-1i*d1(1))*muvec_firstrow(2:M), size(H1,1), size(H1,2)) + ...
        sparse( 1:M-1, M+2:2*M,         -1i*SO_rashba_intrinsic*(d2(2)-1i*d2(1))*muvec_firstrow(1:M-1), size(H1,1), size(H1,2));
        
    
    % expand by spin up/down 
    H0 = [H0 + HSOintrinsic_H0, HSO_rashba_intrinsic_H0; ...
         HSO_rashba_intrinsic_H0b, H0 - HSOintrinsic_H0];
    H1 = [H1 + HSOintrinsic_H1, HSO_rashba_intrinsic_H1; ...
         HSO_rashba_intrinsic_H1b, H1 - HSOintrinsic_H1];
    
end

%% Graphene_Lead_Hamiltonians
%> @brief Creates Hamiltonians H_0 and H_1 of silicene ribbon with zigzag/armchair edge as well as the structure conatining the coordinates of the atomic sites.
%> @image html graphene.jpg
%> @image latex graphene.jpg
%> @param lead_param An instance of structure #Lattice_Graphene (or its subclass) containing the physical parameters.
%> @param M Number of sites in the cross section of the lead.
%> @param End_Type The orientation of the lattice. Set 'A' for lattice with armchair orientation (meaning zigzag edges) or 'Z' for zizgaz orientation (meaning arm-chair edges)
%> @return [1] The Hamiltonian of one slab in the ribbon. 
%> @return [2] The coupling between the slabs. 
%> @return [3] The transverse coupling between the slabs for transverse calculations. 
%> @return [4] A structure Coordinates containing the coordinates of the sites. 
function [H0, H1, H1_transverse, coordinates] = Graphene_Lead_Hamiltonians(obj, lead_param, M, End_Type )
    
    % check the structure containing the physical parameters
    supClasses = superclasses(lead_param);
    if sum( strcmp( supClasses, 'Lattice_Graphene' ) ) == 0
        error(['EQuUs:Lattices:', class(obj), ':Graphene_Lead_Hamiltonians'], 'Invalid type of the input parameter');   
    end
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Graphene_Lead_Hamiltonians'], 'The input parameter M is empty')
    end
    
    epsilon = lead_param.epsilon;
    vargamma = lead_param.vargamma;
    deltaAB = lead_param.deltaAB;
    


    if strcmp(End_Type, 'A')
        [H0,H1, H1_transverse, coordinates] = obj.Armchair_End(epsilon,vargamma,M);
    elseif strcmp(End_Type, 'Z')
        [H0, H1, H1_transverse, coordinates] = obj.Zigzag_End(epsilon,vargamma,M);
    else
        warning('Unrecognized Lead type, creating lead with zigzag edge');
        [H0,H1, H1_transverse, coordinates] = Armchair_End(epsilon,vargamma,M);
    end
    
    % staggered A-B potential tuned by perpendicular electric field   
    H_DeltaAB = sparse( 1:2:M, 1:2:M, deltaAB, size(H0,1), size(H0,2)) + sparse( M+2:2:2*M, M+2:2:2*M, deltaAB, size(H0,1), size(H0,2)) + ...
            sparse( 2:2:M, 2:2:M, -deltaAB, size(H0,1), size(H0,2)) + sparse( M+1:2:2*M, M+1:2:2*M, -deltaAB, size(H0,1), size(H0,2));   
        
    H0 = H0 + H_DeltaAB;
    
%     if ~isempty(q)
%         H0 = H0 + H1_transverse*exp(1i*q) + H1_transverse'*exp(-1i*q);
%     end
    
    
    
    
end


end % methods public

methods ( Access = protected )
    
%> @brief Creates Hamiltonians H_0 and H_1 of silicene ribbon with zigzag edge as well as the structure conatining the coordinates of the atomic sites.
%> @param epsilon The onsite potential.
%> @param vargamma The hopping amplitude.
%> @param M Number of sites in the cross section of the lead.
%> @param q The transverse momentum quantum number.
%> @return [1] The Hamiltonian of one slab in the ribbon. 
%> @return [2] The coupling between the slabs. 
%> @return [3] The transverse coupling between the slabs for transverse calculations. 
%> @return [4] A structure Coordinates containing the coordinates of the sites. 
function [H0,H1,H1_transverse,coordinates] = Armchair_End(obj, epsilon,vargamma,M)
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Armchair_End'], 'The input parameter M is empty')
    end
    
    
    H0 = sparse( 1:2:M-1, 2:2:M, -vargamma, 2*M, 2*M) + ...
        sparse( M+2:2:2*M-1, M+3:2:2*M, -vargamma, 2*M, 2*M) + ...
        sparse( 1:M, M+1:2*M, -vargamma, 2*M, 2*M);
    H0 = H0 + H0' + sparse( 1:2*M, 1:2*M, epsilon, 2*M, 2*M);
    
    H1 = sparse( M+1:2*M, 1:M, -vargamma, 2*M, 2*M);
    
     
     coordinates = structures('coordinates');
     coordinates.a = zeros(2,1);
     coordinates.x = zeros(2*M,1);
     coordinates.y = zeros(2*M,1);
     deltax1 = 2;
     deltax2 = 1;
     deltay  = sqrt(3)/2;
          
    coordinates.x(1:2:M) = ((1:2:M)-1)*(deltax1+deltax2)/2; %in units rCC
     coordinates.x(2:2:M) = ((2:2:M)-2)*(deltax1+deltax2)/2 + deltax2; %in units rCC
     coordinates.x(M+(1:2:M)) = ((1:2:M)-1)*(deltax1+deltax2)/2 - deltax2/2; %in units rCC
     coordinates.x(M+(2:2:M)) = ((2:2:M)-2)*(deltax1+deltax2)/2 + deltax1 - deltax2/2; %in units rCC        
     coordinates.y = [zeros(M,1); ones(M,1)*deltay]; %in units rCC
    coordinates.a = [0; 2*deltay];% 1D lattice constant in uits of rCC 
    
    if mod(M,2) == 0
        H1_transverse = sparse( 2*M, M+1, -vargamma, 2*M, 2*M);
        %H1_transverse = sparse( M+1, 2*M, -vargamma, 2*M, 2*M);
        coordinates.b = [M/2*3; 0];
    else
           
        H1_transverse = [];
    end
        
end

%% Zigzag_End
%> @brief Creates Hamiltonians H_0 and H_1 of silicene ribbon with armchair edge as well as the structure conatining the coordinates of the atomic sites.
%> @param epsilon The onsite potential.
%> @param vargamma The hopping amplitude.
%> @param M Number of sites in the cross section of the lead.
%> @param q The transverse momentum quantum number.
%> @param varargin Optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @return [1] The Hamiltonian of one slab in the ribbon. 
%> @return [2] The coupling between the slabs. 
%> @return [3] The transverse coupling between the slabs for transverse calculations. 
%> @return [4] A structure Coordinates containing the coordinates of the sites. 
function [H0, H1, H1_transverse, coordinates] = Zigzag_End(obj, epsilon,vargamma,M)
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Zigzag_End'], 'The input parameter M is empty')
    end

    coordinates = structures('coordinates');
    coordinates.a = [0;3]; %in units rCC
    
    vargamma_nnn = 0*vargamma;
    H0nnn = sparse(1:M, 2*M+1:3*M, -vargamma_nnn, 4*M, 4*M) + sparse(1:M-1, 2*M+2:3*M, -vargamma_nnn, 4*M, 4*M) + sparse(1:M-1, 2:M, -vargamma_nnn, 4*M, 4*M) +...
            sparse(M+1:2*M, 3*M+1:4*M, -vargamma_nnn, 4*M, 4*M) + sparse(M+2:2*M, 3*M+1:4*M-1, -vargamma_nnn, 4*M, 4*M) + sparse(M+1:2*M-1, M+2:2*M, -vargamma_nnn, 4*M, 4*M) +...
            sparse(2*M+1:3*M-1, 2*M+2:3*M, -vargamma_nnn, 4*M, 4*M) + ...
            sparse(3*M+1:4*M-1, 3*M+2:4*M, -vargamma_nnn, 4*M, 4*M);
        
    H1nnn = sparse(3*M+1:4*M, M+1:2*M, -vargamma_nnn, 4*M, 4*M) + sparse(3*M+1:4*M-1, M+2:2*M, -vargamma_nnn, 4*M, 4*M) + ...
        sparse(2*M+1:3*M, 1:M, -vargamma_nnn, 4*M, 4*M) + sparse(2*M+2:3*M, 1:M-1, -vargamma_nnn, 4*M, 4*M);
    
    H0 = sparse(1:M, M+1:2*M, -vargamma, 4*M, 4*M) + sparse(1:M-1, M+2:2*M, -vargamma, 4*M, 4*M) + ...
        sparse(M+1:2*M, 2*M+1:3*M, -vargamma, 4*M, 4*M) + ...
        sparse(2*M+1:3*M, 3*M+1:4*M, -vargamma, 4*M, 4*M) + sparse(2*M+2:3*M, 3*M+1:4*M-1, -vargamma, 4*M, 4*M);
    H0 = H0 + H0' + sparse( 1:4*M, 1:4*M, epsilon, 4*M, 4*M) + H0nnn + H0nnn';
    
    H1 = sparse(3*M+1:4*M, 1:M, -vargamma, 4*M, 4*M) + H1nnn;
    
    
     coordinates.x = zeros(4*M,1); %in units rCC
     coordinates.y = zeros(4*M,1); %in units rCC
    
    coordinates.x(1:M) = sqrt(3)*(1:M) - sqrt(3)/2;  %in units rCC
    coordinates.y(1:M) = 0;  %in units rCC
    
    coordinates.x(M+1:2*M) = sqrt(3)*(0:M-1);  %in units rCC
    coordinates.y(M+1:2*M) = 0.5;  %in units rCC
    
    coordinates.x(2*M+1:3*M) = sqrt(3)*(0:M-1);  %in units rCC
    coordinates.y(2*M+1:3*M) = 1.5;  %in units rCC
    
    coordinates.x(3*M+1:4*M) = sqrt(3)*(1:M) - sqrt(3)/2;  %in units rCC
    coordinates.y(3*M+1:4*M) = 2;  %in units rCC
    
    %H1_transverse = sparse(M, M+1, -vargamma, 4*M, 4*M) + sparse(4*M, 2*M+1, -vargamma, 4*M, 4*M);
    H1_transverse = sparse(M+1, M, -vargamma, 4*M, 4*M) + sparse(2*M+1, 4*M, -vargamma, 4*M, 4*M);
    coordinates.b = [sqrt(3)*M, 0];
    
end

end % methods protected

end