Siesta_Interface.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - Custom_Hamiltonians
%    Copyright (C) 2015-2018 Peter Rakyta, Ph.D., David Visontai, 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/.
%
%   Interface to import Siesta hamiltonian from Syslabe.HSX provided by other codes or created
%   manually
%%

%> @file Siesta_Interface.m
%> @brief Creates the Hamiltonians and overlap integrals from the loaded data from the Siesta package.
%
%> @Available
%> EQuUs v4.9.... or later
%> @param WorkingDir is declared in the input, this is where the SIESTA files are to be found
%> @param Param contains input parameters
%> @param q k-point, where the Hamiltonian and Overlap matrices should be constructed
%> @return The Leads' H0 and H1 as HL0, SL0, HL1, SL1, the hamiltonian of the scattering part as HSC and SSC and coupling between the scattering region and the leads as Hc and SC. Coordinates and the Fermi energy is returned as well.
function [HL0, SL0, HL1, SL1, coordsLeads, HSC, SSC, coordsSC, ...
    H1_t, S1_t, HSC_t, SSC_t, HC, SC, EF] = Siesta_Interface(WorkingDir, Param, q)

    Syslabel = Param.scatter.FileName(1:end-4);
    XML_file = [WorkingDir, '/', Syslabel, '.xml'];
    HSX_file = [WorkingDir, '/', Syslabel, '.HSX'];
    
    %THESE ARE THE FIELDS FROM XML THAT ARE USED LATER
    fieldsfromxml = {'siesta:E_Fermi', 'siesta:nkpoints', 'Initial System', 'siesta:kpoints', 'siesta:kpt_band'};
    xml_params = Read_Siesta_XML(XML_file, fieldsfromxml, false);
    EF = xml_params.E_Fermi;
    % Use only those k-points that have kz=kpoint(3)=0.0000 ASSUMING THAT
    % THE DIRECTION OF THE TRANSPORT IS ALONG THE Z AXIS
    kpoints = xml_params.kpoints(abs(xml_params.kpoints(:,3))<0.0000001,:);

        % bug fixed: variable atom confronting with object atom
        % see doc. MATLAB/Octave language support for Atom at https://atom.io/packages/language-matlab

    H1_t = [];
    S1_t = [];
    HSC_t = [];
    SSC_t = [];
        
%% CREATE ATOM VECTOR
     NumLeads=length(Param.Leads);
     natoms=length(Param.scatter(1).Atom);
     for ilead=1:NumLeads
          natoms = natoms + length(Param.Leads{ilead}.Atom);
     end

     atom=zeros(natoms,5);

    Number_of_layers_in_leads = zeros(NumLeads);
    Terminating_layer_in_leads = zeros(NumLeads);
    Lead_orientation = zeros(NumLeads);
     for ilead=1:NumLeads
        if size(Param.Leads{ilead}.Atom)~=0;
            Number_of_layers_in_leads(ilead) = max([Param.Leads{ilead}.Atom(1).PL, Param.Leads{ilead}.Atom(end).PL]);
            Terminating_layer_in_leads(ilead) = max([Param.Leads{ilead}.Atom(1).PL, Param.Leads{ilead}.Atom(end).PL]);
            Lead_orientation(ilead) = Param.Leads{ilead}.Lead_Orientation;
        end
          for atm=Param.Leads{ilead}.Atom
               atom(atm.Id,2) = ilead;
               atom(atm.Id,3) = atm.PL;
          end
     end

     atom(:,1)=1:natoms;

%% READ SCATTERING REGION FROM SIESTA HSX

    
    %[HS_EM, coordsSC, H1_t, S1_t, HSC_t, tiorb] = Siesta_Equus_Read_EM(Param, HSX_file, xml_params, q);
    [HS_EM, coordsSC, tiorb, Param, ScatterAtoms] = Read_Siesta_Scatter(Param, HSX_file, xml_params, q);
    
    %THIS SHOWS WHICH ORBITAL BELONGS TO WHICH ATOM
    iorb = [tiorb(:,1), atom(tiorb(:,1),2:3), tiorb(:,2:end)];
    
    % THIS JUST INDICATES WHETHER THE HAMILTONIAN IS IN ONE BLOCK OR TWO
    % IN CASE OF COLLINEAR POLARIZED SPIN CALCULATION IT IS IN 2 SUBBLOCKS 
    % OTHERWISE IN 1 SUBBLOCK
    gspin = 1;
    if Param.scatter.nspin==2
        gspin = 2;
        Param.scatter.ElementNumber = [Param.scatter.ElementNumber, Param.scatter.ElementNumber];
        Param.scatter.ElementSymbol = [Param.scatter.ElementSymbol, Param.scatter.ElementSymbol];
        Param.scatter.Atom = [Param.scatter.Atom, Param.scatter.Atom];
        siorb = size(iorb, 1);
        %iorb = [iorb; iorb];
        %iorb(siorb+1:end, 1) = iorb(1:siorb, 1) + siorb;
    end
% READ SCATTERING INFORMATIONS ---------------------------------------

    
%% ORBITAL INFOS

    %CREATE ORBITAL REGISTER
    Orbitals.Mol = intersect(find(iorb(:,2)==0), find(iorb(:,3)==0));
    Orbitals.ExtMol = Orbitals.Mol';
    for ilead=1:NumLeads
        if Terminating_layer_in_leads(ilead)==1 & ~Param.Leads{ilead}.Use_External_Lead 
            fprintf('There must be more than 1 unit cell for leads if no external is provided!'); 
            quit; 
        end
        for ipl=1:Terminating_layer_in_leads(ilead)
            aux = intersect(find(iorb(:,2)==ilead), find(iorb(:,3)==ipl));
            Orbitals.Lead{ilead}.PL{ipl} = aux; 
            if ipl<Terminating_layer_in_leads(ilead)
                Orbitals.ExtMol = [Orbitals.ExtMol, aux'];
            end
        end
    end
    Orbitals.ExtMol = sort(Orbitals.ExtMol);

% END OF ORBITAL INFOS    ------------------------------------
    
LeadAtoms = {};
%% SET UP LEADS
%    for ispin = 1:gspin
       %for ik = 1:length(param.scatter.kpoints)
         for id_lead=1:NumLeads                        
            if Param.Leads{id_lead}.Use_External_Lead == 1
                %% READ EXTERNAL LEADS IF REQUESTED
                disp(' Reading LEAD Hamiltonian from Siesta')
                               
                %aux = orbLead(ilead,:,leadp(ilead,2));
                aux = Orbitals.Lead{id_lead}.PL{1};
                shift = mean(mean(HS_EM(1).H(aux,aux)));
                [HS_Lead(id_lead).lead, coordsLeads{id_lead}, LeadAtoms{id_lead}] = Read_Siesta_Lead(id_lead, Param, HSX_file, xml_params, shift, q);
                
                HL0{id_lead} = HS_Lead(id_lead).lead(1).H0;
                SL0{id_lead} = HS_Lead(id_lead).lead(1).S0;
                HL1{id_lead} = HS_Lead(id_lead).lead(1).H1;
                SL1{id_lead} = HS_Lead(id_lead).lead(1).S1;
                H1_t{id_lead} = zeros(size(HL1)); 
                S1_t{id_lead} = zeros(size(SL1));

            else
                %% OR FROM SCATTERING REGION
                clear aux0 aux1
                LastTerminatingLayer = Terminating_layer_in_leads(id_lead);
                if LastTerminatingLayer
                    aux0 = Orbitals.Lead{id_lead}.PL{LastTerminatingLayer};
                    aux1 = Orbitals.Lead{id_lead}.PL{LastTerminatingLayer-1};
                    HL0{id_lead} = HS_EM(1).H(aux0, aux0);
                    HL1{id_lead} = HS_EM(1).H(aux0, aux1);
                    SL0{id_lead} = HS_EM(1).S(aux0, aux0);
                    SL1{id_lead} = HS_EM(1).S(aux0, aux1);
                    coordsLeads{id_lead}.x = coordsSC.x([aux0; aux1]);
                    coordsLeads{id_lead}.y = coordsSC.y([aux0; aux1]);
                    coordsLeads{id_lead}.z = coordsSC.z([aux0; aux1]);
                    LeadAtoms{id_lead} = ScatterAtoms([aux0; aux1]);

                    if Param.scatter.nspin==2
                        HL0 = [HL0{id_lead}, sparse([],[],[], size(HL0{id_lead},1), size(HL0{id_lead},2)); sparse([],[],[], size(HL0{id_lead},1), size(HL0{id_lead},2)), HS_EM(2).H(aux0, aux0)];
                        HL1 = [HL1{id_lead}, sparse([],[],[], size(HL1{id_lead},1), size(HL1{id_lead},2)); sparse([],[],[], size(HL1{id_lead},1), size(HL1{id_lead},2)), HS_EM(2).H(aux0, aux1)];
                        coordsLeads{id_lead}.x = [coordsLeads{id_lead}.x, coordsLeads{id_lead}.x];
                        coordsLeads{id_lead}.y = [coordsLeads{id_lead}.y, coordsLeads{id_lead}.y];
                        coordsLeads{id_lead}.z = [coordsLeads{id_lead}.z, coordsLeads{id_lead}.z];
                        LeadAtoms{id_lead} = [LeadAtoms{id_lead}, LeadAtoms{id_lead}];
                        coordsLeads{id_lead}.spinup = [ true(size(HL0{id_lead},1),1); false(size(HL0{id_lead},1),1)];
                    end
%                 else
%                     HL0{id_lead} = [];
%                     HL1{id_lead} = [];
%                     SL0{id_lead} = [];
%                     SL1{id_lead} = [];
%                     coordsLeads{id_lead}.x = [];
%                     coordsLeads{id_lead}.y = [];
%                     coordsLeads{id_lead}.z = [];
                end

                    H1_t{id_lead} = zeros(size(HL1)); 
                    S1_t{id_lead} = zeros(size(SL1));
          end
         end; 
%    end; 
    
% END OF LEADS SETUP INFOS    ------------------------------------    
    


    coordsSC.x = coordsSC.x(Orbitals.Mol);
    coordsSC.y = coordsSC.y(Orbitals.Mol);
    coordsSC.z = coordsSC.z(Orbitals.Mol);
    ScatterAtoms = ScatterAtoms(Orbitals.Mol);

    if Param.scatter.nspin==2
        coordsSC{id_lead}.x = [coordsSC{id_lead}.x, coordsSC{id_lead}.x];
        coordsSC{id_lead}.y = [coordsSC{id_lead}.y, coordsSC{id_lead}.y];
        coordsSC{id_lead}.z = [coordsSC{id_lead}.z, coordsSC{id_lead}.z];
        coordsSC{id_lead}.Atoms = [coordsSC{id_lead}.Atoms, coordsSC{id_lead}.Atoms];
        coordsSC{id_lead}.spinup = [ true(size(HL0{id_lead},1),1); false(size(HL0{id_lead},1),1)];
    end
    
    % FROM THE WHOLE HAMILTONIAN CUT OUT THE EFFECT OF PERIODICITY
 for ispin = 1:gspin
    for il1=1:NumLeads
        for il2=il1+1:NumLeads
            Leadorb1 = [];
            for ipl1=1:Terminating_layer_in_leads(il1)
               %Leadorb1 = [Leadorb1; Orbitals.Lead{il1}.PL{ipl1}];
                Leadorb1 = [Leadorb1, Orbitals.Lead{il1}.PL{ipl1}];
                Leadorb2 = [];
                for ipl2=1:Terminating_layer_in_leads(il2)
                    %Leadorb2 = [Leadorb2; Orbitals.Lead{il2}.PL{ipl2}];
                    Leadorb2 = [Leadorb2, Orbitals.Lead{il2}.PL{ipl2}'];
                    HS_EM(ispin).H(Leadorb1, Leadorb2) = 0;
                    HS_EM(ispin).H(Leadorb2, Leadorb1) = 0;
                end
            end
        end
    end
  end
    
    

        %% LOAD SCATTERING REGION INTO HSC, SSC    
        %ez jo egyebkent
        HSC = HS_EM(1).H(Orbitals.Mol, Orbitals.Mol);
        SSC = HS_EM(1).S(Orbitals.Mol, Orbitals.Mol);
          if Param.scatter.nspin==2
            HSC = [HSC, sparse([],[],[], size(HSC, 1), size(HSC, 2)); sparse([],[],[], size(HSC, 1), size(HSC, 2)), HS_EM(2).H(Orbitals.Mol, Orbitals.Mol)];
            SSC = [SSC, sparse([],[],[], size(SSC, 1), size(SSC, 2)); sparse([],[],[], size(SSC, 1), size(HSC, 2)), HS_EM(2).S(Orbitals.Mol, Orbitals.Mol)];
          end
        %% LOAD INTERFACE INTO HC, SC
        for ilead =1:NumLeads
            LastTerminatingLayer = Terminating_layer_in_leads(ilead);
            if LastTerminatingLayer 
                C = Orbitals.Lead{ilead}.PL{max(LastTerminatingLayer-1, 1)};
                %ez jo!!!
                HC{ilead} = HS_EM(ispin).H(C, Orbitals.Mol);
                SC{ilead} = HS_EM(ispin).S(C, Orbitals.Mol);
            else
                HC{ilead} = [];
                SC{ilead} = [];
            end
        end

   HSC_t = zeros(size(HSC));
   SSC_t = zeros(size(SSC));
   
end %function




    
%     %Returns with element number
%     function indx = get_atomindex(sym)
%     Syms={'H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si','P','S','Cl',...
%         'Ar','K','Ca','Sc','Ti','V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',...
%         'Se','Br','Kr','Rb','Sr','Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In',...
%         'Sn','Sb','Te','I','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb',...
%         'Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb',...
%         'Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa','U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No',...
%         'Lr','Rf','Db','Sg','Bh','Hs','Mt','Ds','Rg','Cn','Uut','Uuq','Uup','Uuh','Uus','Uuo'};
% 
%         indx=find(ismember(Syms,sym)==1);
% 
% end