Graphene_Bilayer_Lead_Hamiltonians.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - Graphene_Bilayer_Lead_Hamiltonians
%    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 lattices Lattices
%> @{
%> @file Graphene_Bilayer_Lead_Hamiltonians.m
%> @brief A class to create the Hamiltonian of one unit cell in a translational invariant graphene bilayer lead.
%> @} 
%> @brief A class to create the Hamiltonian of one unit cell in a translational invariant graphene bilayer lead.
%%
classdef Graphene_Bilayer_Lead_Hamiltonians < Hex_Lead_Hamiltonians
    

methods ( Access = public )
    
 

%% Graphene_Bilayer_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_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_Hamiltonians'], 'Invalid type of the input parameter');   
    end
    
    if isempty(M)
        error(['EQuUs:Lattices:', class(obj), ':Graphene_Bilayer_Hamiltonians'], 'The input parameter M is empty')
    end    
    
    [H0_lower,H1_lower,H1_transverse_lower,coordinates_lower] = obj.Graphene_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_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
    
    % combining sites of the upper and lower layer
    coordinates = coordinates_lower.Combine( coordinates_upper );
    
    
    
end



%% Graphene_Bilayer_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_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_Hamiltonians2'], 'Invalid type of the input parameter');   
    end
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Graphene_Bilayer_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_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);
        % Sites with false values will be removed from the structure.
        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
    
        % removing sites from the coordinates
        coordinates = coordinates.KeepSites( logical_indexes );
        
    
    
    % armchair edged bilayer ribbon
    elseif strcmp(End_Type, 'Z')
        [H0, H1, H1_transverse, coordinates] = obj.Graphene_Bilayer_Hamiltonians(lead_param, M, End_Type, varargin{:});        
    end
    
    
end



%% Graphene_Bilayer_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_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_Hamiltonians3'], 'Invalid type of the input parameter');   
    end
    
    if isempty(M)
        error(['EQuUs:', class(obj), ':Graphene_Bilayer_Hamiltonians3'], 'The input parameter M is empty')
    end
        
    Mtmp = M+2;
    
    [H0, H1, H1_transverse, coordinates] = obj.Graphene_Bilayer_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);
    % Sites with false values will be removed from the structure.
    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
    
    % removing sites from the coordinates
    coordinates = coordinates.KeepSites( logical_indexes );
    
    
    
end



end % methods public


end