test_graphene_kpoints_siesta_SZ.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_Graphene_kpoint_siesta.m
%> @brief Testfile to check SIESTA (dft) interface on a graphene strip with k-point sampling.
%> @Available
%> EQuUs v4.9.... or later
%> @}
function test_graphene_kpoints_siesta_SZ( filenum )

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

    Label = 'graphene_kpoints_siesta_SZ';

    % create input structures
    filenameIn  = ['Input_',Label,'.xml'];
    filenameOut = ['output-',Label,'.xml'];
    [Opt, param] = parseInput( filenameIn );

    Silent=1;

    %THESE ARE THE FIELDS FROM XML THAT ARE USED LATER
    fieldsfromxml = {'siesta:kpoints', 'siesta:E_Fermi', 'siesta:nkpoints', 'Initial System', 'siesta:kpoints', 'siesta:kpt_band'};
    siestaxml = [Opt.WorkingDir, '/', param.scatter.FileName(1:end-4),'.xml'];
    xml_params = Read_Siesta_XML(siestaxml, fieldsfromxml);
    kpoints = xml_params.kpoints(:,1:3);
    kweights = xml_params.kpoints(kpoints(:,3)<1e-9, 4);  %only kpoints that are zero in the direction of the transport
    kweights = kweights/sum(kweights);
     Bohr   =  0.529177;                       % Conversion from Bohr to Angstrom

    kpoints = kpoints(kpoints(:,3)<1e-9,:)/Bohr;  % For some reason Siesta stores K-point everywhere in 1/Bohr
    
    Evec = [-3.2];
    check_conductivity = [1.533333333333339 ];
    
    format long;     
    
    disp( '***** test of SIESTA interface, Transport of graphene system with multiple transversal K-points *****' );    

    Conductivity = NaN(length(Evec),1);
    Open_channels1 = NaN(length(Evec),1);
    Open_channels2 = NaN(length(Evec),1);
    
    for idx = 1:length(Evec)
        [Conductivity(idx), Open_channels1(idx), Open_channels2(idx)] = CalculateTransport(Evec(idx));
    end

    function [T, ny1, ny2] = CalculateTransport(E)     
        
        num = length(kpoints);
        T_k = NaN(num,1); 
        ny1_k = NaN(num,1); 
        ny2_k = NaN(num,1); 
        for ikp=1:length(kpoints)
            kpoint = kpoints(ikp, 1:3);
            kweight = kweights(ikp);

            cTwoTerminal = NTerminal('filenameIn', fullfile( pwd, filenameIn), 'filenameOut', fullfile( filenameOut), 'Silent',1, 'q', kpoint);

            DeltaC_tmp = 0;
            [Conductivity_tmp, ny_tmp, DeltaC_tmp] = cTwoTerminal.Transport(E, 'constant_channels', false, 'SelfEnergy', 1);
      
            T_k(ikp) = -Conductivity_tmp(1,1)*kweight;
            ny1_k(ikp) = ny_tmp(1)*kweight;
            ny2_k(ikp) = ny_tmp(2)*kweight;
            
        end      
        T = sum(T_k);
        ny1 = sum(ny1_k);
        ny2 = sum(ny2_k);
    end
    
    % checking conductances
    checkpoint = max(Conductivity - check_conductivity);
    if checkpoint > 1e-7
        warning(['EQuUs:Tests:', fncname, ':checkpoint failed with error ', num2str( checkpoint)]);
        Conductivity
    end


end