Read_Siesta_Scatter.m

Go to the documentation of this file.

%> @file Read_Siesta_Scatter.m
%> @brief Creates the Hamiltonians and Overlap matrices from the HSX file data of the SIESTA output.
%> @Available
%> EQuUs v4.9 or later
% 
%> @brief Creates the Hamiltonians and Overlap matrices from the HSX file data of the SIESTA output.
%> @param Param contains general input parameters
%> @param SIESTA_HSX_file is the file, where the hamiltonian is stored in a unique sparse format of SIESTA
%> @param SIESTA_XML_params parameters from SIESTA xml output
%> @param q k-point, where the Hamiltonian and Overlap matrices should be constructed
%> @return The constructed Hamiltonian of the whole system as HS, atomic coordinates and a description of the orbital system
    function [HS, coordsSC, iorb, Param, Atoms] = Read_Siesta_Scatter(Param, SIESTA_HSX_file, SIESTA_XML_params, q)

     %SOME CONSTANTS USE IN SIESTA
     RytoeV = 1/13.60580;                      % in SIESTA_HSX_file values are stored in Ry and Bohr
     Bohr   =  0.529177;                       % These are conversion units to eV and Angstrom

    %% TODO kpoints needed?
    %param.scatter.kpoints = xml_params.kpoints;
    Param.scatter.ElementNumber = SIESTA_XML_params.elmtxv;
    Param.scatter.ElementSymbol = SIESTA_XML_params.spxv;
    Param.scatter.LatticeVectors = SIESTA_XML_params.Lattice_vectors;
        
     %COORDS and %LATTICE VECTORS
    coordsSC = structures('coordinates');
    
    % TEMP
    HSXmat= 'SIESTAHSX.mat';
    
     %% READ FILE "system_label.HSX"
    if exist(HSXmat, 'file') == 0 | true
        hsx = read_hsx(SIESTA_HSX_file);
        
        %IF USING THE FORTRAN HSX READER, THEN VARIABLES HAVE TO BE CONVERTED TO DOUBLES
        hsx.numh=double(hsx.numh);
        hsx.listhptr=double(hsx.listhptr);
        hsx.indxuo=double(hsx.indxuo);

        no_u     = double(hsx.no_u);
        nspin    = double(hsx.nspin);
        Param.scatter.nspin = hsx.nspin;

        %Set units to eV
        hsx.hamilt = hsx.hamilt/RytoeV;
        
        Atoms = zeros(hsx.no_u);

        %% CREATE ORBITAL INFO
        iorb = [];
        for iuo=1:hsx.no_u
            iat    = hsx.iaorb(iuo);
            iphorb = hsx.iphorb(iuo);
            spec   = hsx.isa(iat);
            element_index = Param.scatter.ElementNumber(iat);
            coordsSC.x(iuo) = SIESTA_XML_params.coordinates(iat,1);
            coordsSC.y(iuo) = SIESTA_XML_params.coordinates(iat,2);
            coordsSC.z(iuo) = SIESTA_XML_params.coordinates(iat,3);
            Atoms(iuo) = iat;
            dum    = [ iat, hsx.nquant(spec,iphorb), hsx.lquant(spec,iphorb), 0, hsx.zeta(spec,iphorb),...
                element_index, ...
                Get_Atomic_Mass(element_index)];
            iorb = [iorb; dum];
        end

        %% IF IT IS FROM A COLLINEAR SPIN POLARIZED CALCULATION THEN 
        % WE DOUBLE THE SIZE OF THE HAMILTONIAN AND STORE THEM IN SEPARATE BLOCKS
        % BUT IN A NONCOLLINEAR CASE 
        if nspin > 2
            iorb = [iorb; iorb];
            coordsSC.x = [coordsSC.x; coordsSC.x];
            coordsSC.y = [coordsSC.y; coordsSC.y];
            coordsSC.z = [coordsSC.z; coordsSC.z];
            element_index = [element_index; element_index];
            Atoms = [Atoms; Atoms];
        end

        %% Find supercell size
        % Transport direction is always into the direction of the z-coordinate
        sux = floor( max(hsx.xij(1,:)) / Param.scatter.LatticeVectors(1,1)+0.5);
        sux = 2 * sux + 1;
        suy = floor( max(hsx.xij(2,:)) / Param.scatter.LatticeVectors(2,2)+0.5);
        suy = 2 * suy + 1;
        suz = floor( max(hsx.xij(3,:)) / Param.scatter.LatticeVectors(3,3)+0.5);
        if suz>0
            throw('It should not have kpoints into z directions, which is the transport direction')
        end

        %COMPOSING THE Hk FROM HAMILTONIAN ELEMENTS
        %EM PART, THE SCATTERING REGION
        %IN CASE OF JOSEPHSON CURRENT MODE, THIS EXCLUDES THE INTERFACE
        %WHICH IS THE PART OF THE ELECTRODES, CLOSEST TO EM, THAT IS NOT USED 
        %FOR THE LEADS CALCULATION
        tt=[]; tk=[];
        for iuo=1:no_u
            tt = [ tt 1:hsx.numh(iuo)];
            tk= [tk hsx.listhptr(iuo).*ones(hsx.numh(iuo),1)'];%'
        end
        ind=(tt+tk)';%'

        iuo=[];
        for ii=1:no_u
            iuo = [iuo double(ii).*ones(hsx.numh(ii),1)'];%'
        end
        jo=hsx.listh(ind);
        juo=(hsx.indxuo(jo))';%'

        % Convert hsx.xij from Bohr to Angstrom
        xij = hsx.xij*Bohr;

        %% TODO maybe coordinates(iuo',1) = coordsSC.x(iorb(iuo',1)
        Rbloch(:,1) = xij(1,ind)' + coordsSC.x(iuo)' - coordsSC.x(juo)';
        Rbloch(:,2) = xij(2,ind)' + coordsSC.y(iuo)' - coordsSC.y(juo)';
        Rbloch(:,3) = xij(3,ind)' + coordsSC.z(iuo)' - coordsSC.z(juo)';
        %xij = hsx.xij;
        %save(HSXmat, 'hsx', 'coordsSC', 'iorb', 'Rbloch', 'ind', 'iuo', 'juo', 'no_u');
    else
        load(HSXmat);
    end

    %% Find supercell size (experiment with this later)
    %sux = floor( max(hsx.xij(1,:)) / Param.scatter.LatticeVectors(1,1)+0.5);
    %sux = 2 * sux + 1;
    %suy = floor( max(hsx.xij(2,:)) / Param.scatter.LatticeVectors(2,2)+0.5);
    %suy = 2 * suy + 1;
    %suz = floor( max(hsx.xij(3,:)) / Param.scatter.LatticeVectors(3,3)+0.5);
    %if suz>0
    %    throw('It should not have kpoints into z directions, which is the transport direction')
    %end
    
     %COMPOSING THE Hk FROM HAMILTONIAN ELEMENTS
     %EM PART, THE SCATTERING REGION
     %IN CASE OF JOSEPHSON CURRENT MODE, THIS EXCLUDES THE INTERFACE
     %WHICH IS THE PART OF THE ELECTRODES, CLOSEST TO EM, THAT IS NOT USED 
     %FOR THE LEADS CALCULATION
    
     HS=[];
    if length(q)>1
        eikr = exp(-1.0i* (q(1).*Rbloch(ind,1)+q(2).*Rbloch(ind,2)+q(3).*Rbloch(ind,3)));
    else
        eikr = ones(length(ind), 1);
    end
    if Param.scatter.nspin <= 2
        HS(1).S = sparse(juo,iuo, hsx.Sover.*eikr,no_u,no_u);
        for is=1:Param.scatter.nspin
            HS(is).H = sparse(juo,iuo, hsx.hamilt(:,is).*eikr,no_u,no_u);
        end
    elseif Param.scatter.nspin == 4
        HS(1).S = kron(eye(2),sparse(juo,iuo, hsx.Sover.*eikr,no_u,no_u));
        H1=kron([1  0; 0 0], sparse(juo,iuo, hsx.hamilt(:,1).*eikr,no_u,no_u));
        H2=kron([0  0; 0 1], sparse(juo,iuo, hsx.hamilt(:,2).*eikr,no_u,no_u));
        H3=kron([0  1; 0 0], sparse(juo,iuo, (hsx.hamilt(:,3)-1i*hsx.hamilt(:,4)).*eikr,no_u,no_u));
        H4=kron([0  0; 1 0], sparse(juo,iuo, (hsx.hamilt(:,3)+1i*hsx.hamilt(:,4)).*eikr,no_u,no_u));
        HS(1).H = H1 + H2 + H3 + H4; 
    elseif Param.scatter.nspin == 8
        HS(1).S = kron(eye(2),sparse(juo,iuo, hsx.Sover.*eikr,no_u,no_u));
        H1=kron([1  0; 0 0], sparse(juo,iuo, (hsx.hamilt(:,1)+1i*hsx.hamilt(:,5)).*eikr,no_u,no_u));
        H2=kron([0  0; 0 1], sparse(juo,iuo, (hsx.hamilt(:,2)+1i*hsx.hamilt(:,6)).*eikr,no_u,no_u));
        H3=kron([0  1; 0 0], sparse(juo,iuo, (hsx.hamilt(:,3)-1i*hsx.hamilt(:,4)).*eikr,no_u,no_u));
        H4=kron([0  0; 1 0], sparse(juo,iuo, (hsx.hamilt(:,7)+1i*hsx.hamilt(:,8)).*eikr,no_u,no_u));
        HS(1).H = H1 + H2 + H3 + H4; 
    end
     %end

end %function