test_LatticeTriangle.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_LatticeTriangle.m
%> @brief Testfile to check functionalities related to calculations performed on triangle lattice
%> @Available
%> EQuUs v4.9 or later
%> @}
%> @brief Testfile to check functionalities related to calculations performed on triangle lattice
function test_LatticeTriangle()
    
filename = mfilename('fullpath');
[directory, fncname] = fileparts( filename );

    % creating the shape of the scattering region
    shape=structures('shape');
    shape.width=1;
    shape.height=10;

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

    % creating class CreateLeadHamiltonians
    cCreateLeadHamiltonians = CreateLeadHamiltonians(Opt, param, 'Hanyadik_Lead',[]);
    
    % creating Hamiltonians
    cCreateLeadHamiltonians.CreateHamiltonians();
    %cCreateLeadHamiltonians.CalcSpektrum('ka_min', -1, 'ka_max', 1, 'ka_num', 100, 'toPlot', true)

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

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

    % special points
    GammaPoint = [0;0];
    KPoint_p = (2*b1-b2)/3;
    KPoint_m = -(2*b1-b2)/3;
    MPoint = b1/2;

    %% test of the Fourier transformed Hamiltonians
    disp( '***** test the energy eigenvalues at specific k-points *****' );

    % Calculate the energy eigenvalue at a specific K-point
    kvector = GammaPoint;
    % using EQuUs
    Spectrum_EQuUs = cCreateLeadHamiltonians.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);

    % using the Fourier-transformed Hamiltonian
    Spectrum_Fourier = SpectrumCalc( kvector );
    
    checkpoint = abs( Spectrum_EQuUs - Spectrum_Fourier );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
    % Calculate the energy eigenvalue at a specific K-point
    kvector = KPoint_p;
    % using EQuUs
    Spectrum_EQuUs = cCreateLeadHamiltonians.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);

    % using the Fourier-transformed Hamiltonian
    Spectrum_Fourier = SpectrumCalc( kvector );
    
    checkpoint = abs( Spectrum_EQuUs - Spectrum_Fourier );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
    
    % Calculate the energy eigenvalue at a specific K-point
    kvector = KPoint_m;
    % using EQuUs
    Spectrum_EQuUs = cCreateLeadHamiltonians.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);

    % using the Fourier-transformed Hamiltonian
    Spectrum_Fourier = SpectrumCalc( kvector );
    
    checkpoint = abs( Spectrum_EQuUs - Spectrum_Fourier );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end
    
    
    % Calculate the energy eigenvalue at a specific K-point
    kvector = MPoint;
    % using EQuUs
    Spectrum_EQuUs = cCreateLeadHamiltonians.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);

    % using the Fourier-transformed Hamiltonian
    Spectrum_Fourier = SpectrumCalc( kvector );
    
    checkpoint = abs( Spectrum_EQuUs - Spectrum_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 *****' );
    
    % Calculate the energy eigenvalue at a specific K-point
    kvector = 0.3*KPoint_p + 0.7*MPoint;
    
    
    % setting the width of the first lead
    param.Leads{1}.M = 3;

    % creating class CreateLeadHamiltonians
    cCreateLeadHamiltonians_Lead = CreateLeadHamiltonians(Opt, param, 'Hanyadik_Lead',1);
    
    % creating Hamiltonians
    cCreateLeadHamiltonians_Lead.CreateHamiltonians();
    
    % Creating the Hamiltonian of the M width lead.
    Hamiltonian_M3 = cCreateLeadHamiltonians_Lead.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2*param.Leads{1}.M);
    
    % Creating the Hamiltonian of the M=1 width lead.
    Hamiltonian_M1 = cCreateLeadHamiltonians.MomentumDependentHamiltonian( kvector'*a1, kvector'*a2);
    
    difference = eig(full(Hamiltonian_M3)) - ones(param.Leads{1}.M,1)*eig(full(Hamiltonian_M1));
    difference = sort( abs(difference), 'ascend' );
    
    checkpoint = abs( difference(1) );
    if checkpoint > 1e-10
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
    end    
    

%% SpectrumCalc
%> @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 = SpectrumCalc( kvector )
        parameters = param.scatter;
        epsilon = parameters.epsilon;
        vargamma = parameters.vargamma;
        
        ret = epsilon + 2*vargamma*( cos(kvector'*delta1) + cos(kvector'*delta2) + cos(kvector'*delta3) );        
        
    end

end