Read_Siesta_Lead.m

Go to the documentation of this file.

%> @file Read_Siesta_Lead.m
%> @brief Creates the Hamiltonians and overlap integrals from the HSX file data of the Siesta package for the lead
%> @Available
%> EQuUs v4.9 or later
%> @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 shift, the coordinates to fit to the scattering region
%> @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, coordsLeads, Atoms] = Siesta_Equus_Read_Lead(id_lead, Param, SIESTA_HSX_file, SIESTA_XML_params,shift, q)

    %SOME CONSTANTS USE IN SIESTA
     RytoeV = 1/13.60580;                      % 
     Bohr   =  0.529177;                       % Angstrom
    
    Syslabel = Param.Leads{id_lead}.FileName(1:end-4);
    xmlfile = [Syslabel,'.xml'];
    
     %THESE ARE THE FIELD FROM XML THAT ARE USED LATER
     fieldsfromxml = { 'siesta:E_Fermi', 'siesta:nkpoints', 'Initial System' ,'siesta:kpoints','siesta:kpt_band', 'siesta:kscell'};
     Lead_params = Read_Siesta_XML(xmlfile, fieldsfromxml, true);
    Param.Leads{id_lead}.E_Fermi = Lead_params.E_Fermi;
    
    %% EZ ITT RONDA, DE EZ VAN. EZ AZ ADAT NINCS MEG MASHOL
    %% FIND THE DIRECTORY WHERE THE "filename.out" FILE IS PLACED AND READ IT
    output_file =  [Syslabel,'.out'];
%    directory = inputP.Path_Leads{lead}(1:end-cellfun(@length, output_file(end)));
    logfile=fileread(output_file);

    %% FIND SYSTEM LABEL
    foundat = findstr(logfile,'SystemLabel');
    dummy=strsplit(logfile(foundat:foundat+100));
    system_label = dummy(2);
    
    %% FIND SUPERCELL SIZE
    foundat = findstr(logfile,'supercell');
    if ~ isempty(foundat)
            dummy = strsplit(logfile(foundat:foundat+100));
            if ~isempty(str2num(dummy{3})) && ~isempty(str2num(dummy{4})) && ~isempty(str2num(dummy{5}))
             nsc = [str2num(dummy{3}) str2num(dummy{4}) str2num(dummy{5})];
            elseif ~isempty(str2num(dummy{2})) && ~isempty(str2num(dummy{4})) && ~isempty(str2num(dummy{6}))
             nsc = [str2num(dummy{2}) str2num(dummy{4}) str2num(dummy{6})];
            else
             warning(' ************* OH NO! SIESTA OUTPUT CHANGED AGAIN?! ************')
             exit
            end
    else
        nsc = [1 1 1];
    end
    %% SUPERCELL SIZE
    nnsc=nsc(1)*nsc(2);

    
    
    %% 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;
    
    kpoints = Lead_params.kpoints;
    elmtxv = Lead_params.elmtxv;
    spxv = Lead_params.spxv;
    LatticeVectors = Lead_params.Lattice_vectors;
    
    %STRUCTURES FOR COORDINATES
    coordsLeads = structures('coordinates');

    %% SELECT ONLY TRANSVERSE K-POINTS & RESTATE WEIGHTS ??
    %obj.param.Leads{}.kpoints;

     %% READ FILE "system_label.HSX"
     hsx = read_hsx( Param.Leads{id_lead}.FileName );

     %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);
     gspin = hsx.nspin;
     if nspin > 2
       gspin = 1;
     end

     %Set units to eV
     hsx.hamilt = hsx.hamilt/RytoeV;
    
     %COORDS and %LATTICE VECTORS
    coordsLead = structures('coordinates');
    Atoms = zeros(no_u);
    
     %% CREATE ORBITAL INFO
     iorb = [];
    coords = zeros(no_u,3);
     for iuo=1:no_u
          iat    = hsx.iaorb(iuo);
          iphorb = hsx.iphorb(iuo);
          spec   = hsx.isa(iat);
%            element_index = Param.scatter.ElementNumber(iat);
          element_index  = Lead_params.elmtxv(iat);
          coordsLead.x(iuo) = Lead_params.coordinates(iat,1);
          coordsLead.y(iuo) = Lead_params.coordinates(iat,2);
          coordsLead.z(iuo) = Lead_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
    Param.Leads{id_lead}.coordinates = coordsLead;
     %IF IT IS COLLINEAR SPIN POLARIZED THEN WE DOUBLE THE SIZE OF THE
     %HAMILTONIAN AND STORE THEM IN SEPARATE BLOCKS
     if nspin > 2
          iorb = [iorb; iorb];
    end
     

     %COMPOSING THE Hk FROM HAMILTONIAN ELEMENTS
     %EM PART, THE SCATTERING REGION
     %IN CASE OF JOSEPHSON CURRENT MODE, THIS EXCLUDES THE INTERFACE
     %WHICH IS PART OF THE ELECTRODES, CLOSEST TO EM, THAT IS NOT USED 
     %FOR THE LEADS CALCULATION

     %% FIND INDICES BELONGING TO H0 OR H1
     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))';%'

     ind0 = find(jo <= no_u*nnsc);
     jo0  = jo(ind0);
     juo0 = (hsx.indxuo(jo0))';%'
     iuo0 = iuo(ind0);
     
     ind1 = intersect(intersect(find(jo > no_u*nnsc),find(jo <= 2*no_u*nnsc)),find(hsx.xij(3,ind) >= 0.0));
     jo1  = jo(ind1);
     juo1 = (hsx.indxuo(jo1))';%'
     iuo1 = iuo(ind1);

    % Convert hsx.xij from Bohr to Angstrom
    xij = hsx.xij*Bohr;
    
     %% COMPUTE H and S FOR GIVEN LEAD
    Rbloch(:,1) = xij(1,ind)' + coordsLead.x(iuo)' - coordsLead.x(juo)';
    Rbloch(:,2) = xij(2,ind)' + coordsLead.y(iuo)' - coordsLead.y(juo)';
    Rbloch(:,3) = xij(3,ind)' + coordsLead.z(iuo)' - coordsLead.z(juo)';

     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).S0 = sparse(juo0,iuo0, hsx.Sover(ind0).*eikr(ind0),no_u,no_u);
        HS(1).S1 = sparse(juo1,iuo1, hsx.Sover(ind1).*eikr(ind1),no_u,no_u);
        for is=1:Param.scatter.nspin
            HS(is).H0 = sparse(juo0,iuo0, hsx.hamilt(ind0,is).*eikr(ind0),no_u,no_u);
            HS(is).H1 = sparse(juo1,iuo1, hsx.hamilt(ind1,is).*eikr(ind1),no_u,no_u);
        end
    elseif Param.scatter.nspin == 4
        HS(1).S0 = kron(eye(2),sparse(juo0,iuo0, hsx.Sover(ind0).*eikr(ind0),no_u,no_u));
        HS(1).S1 = kron(eye(2),sparse(juo1,iuo1, hsx.Sover(ind1).*eikr(ind1),no_u,no_u));
        H1=kron([1  0; 0 0], sparse(juo0,iuo0, hsx.hamilt(ind0,1).*eikr(ind0),no_u,no_u));
        H2=kron([0  0; 0 1], sparse(juo0,iuo0, hsx.hamilt(ind0,2).*eikr(ind0),no_u,no_u));
        H3=kron([0  1; 0 0], sparse(juo0,iuo0, (hsx.hamilt(ind0,3)-1i*hsx.hamilt(ind0,4)).*eikr(ind0),no_u,no_u));
        H4=kron([0  0; 1 0], sparse(juo0,iuo0, (hsx.hamilt(ind0,3)+1i*hsx.hamilt(ind0,4)).*eikr(ind0),no_u,no_u));
        HS(1).H0 = H1 + H2 + H3 + H4; 
        H1=kron([1  0; 0 0], sparse(juo1,iuo1, hsx.hamilt(ind1,1).*eikr(ind1),no_u,no_u));
        H2=kron([0  0; 0 1], sparse(juo1,iuo1, hsx.hamilt(ind1,2).*eikr(ind1),no_u,no_u));
        H3=kron([0  1; 0 0], sparse(juo1,iuo1, (hsx.hamilt(ind1,3)-1i*hsx.hamilt(ind1,4)).*eikr(ind1),no_u,no_u));
        H4=kron([0  0; 1 0], sparse(juo1,iuo1, (hsx.hamilt(ind1,3)+1i*hsx.hamilt(ind1,4)).*eikr(ind1),no_u,no_u));
        HS(1).H1 = H1 + H2 + H3 + H4; 
    elseif Param.scatter.nspin == 8
        HS(1).S0 = kron(eye(2),sparse(juo0,iuo0, hsx.Sover(ind0).*eikr(ind0),no_u,no_u));
        HS(1).S1 = kron(eye(2),sparse(juo1,iuo1, hsx.Sover(ind1).*eikr(ind1),no_u,no_u));
        H1=kron([1  0; 0 0], sparse(juo0,iuo0, (hsx.hamilt(ind0,1)+1i*hsx.hamilt(ind0,5)).*eikr(ind0),no_u,no_u));
        H2=kron([0  0; 0 1], sparse(juo0,iuo0, (hsx.hamilt(ind0,2)+1i*hsx.hamilt(ind0,6)).*eikr(ind0),no_u,no_u));
        H3=kron([0  1; 0 0], sparse(juo0,iuo0, (hsx.hamilt(ind0,3)-1i*hsx.hamilt(ind0,4)).*eikr(ind0),no_u,no_u));
        H4=kron([0  0; 1 0], sparse(juo0,iuo0, (hsx.hamilt(ind0,7)+1i*hsx.hamilt(ind0,8)).*eikr(ind0),no_u,no_u));
        HS(1).H0 = H1 + H2 + H3 + H4; 
        H1=kron([1  0; 0 0], sparse(juo1,iuo1, (hsx.hamilt(ind1,1)+1i*hsx.hamilt(ind1,5)).*eikr(ind1),no_u,no_u));
        H2=kron([0  0; 0 1], sparse(juo1,iuo1, (hsx.hamilt(ind1,2)+1i*hsx.hamilt(ind1,6)).*eikr(ind1),no_u,no_u));
        H3=kron([0  1; 0 0], sparse(juo1,iuo1, (hsx.hamilt(ind1,3)-1i*hsx.hamilt(ind1,4)).*eikr(ind1),no_u,no_u));
        H4=kron([0  0; 1 0], sparse(juo1,iuo1, (hsx.hamilt(ind1,7)+1i*hsx.hamilt(ind1,8)).*eikr(ind1),no_u,no_u));
        HS(1).H1 = H1 + H2 + H3 + H4; 
    end

%% Change Lead orientation if needed
     if Param.Leads{id_lead}.Lead_Orientation == -1;
        HS(1).S1 = HS(1).S1';
        for is=1:gspin
           HS(is).H1 = HS(is).H1';
         end
    end

    %% Correction
     correction = mean(mean(HS(1).H0))-shift;

       %for is=1:obj.param.scatter.nspin
      for is=1:gspin
         HS(is).H0 = HS(is).H0 - correction*HS(1).S0;
         HS(is).H1 = HS(is).H1 - correction*HS(1).S1;
      end
      

end %function