Custom_Hamiltonians.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - Custom_Hamiltonians
%    Copyright (C) 2015-2016 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/.
%
%> @addtogroup basic Basic Functionalities
%> @{
%> @file Custom_Hamiltonians.m
%> @brief A class to import custom Hamiltonians provided by other codes or created manually
%> @} 
%> @brief A class to import custom Hamiltonians provided by other codes or created manually
%> @Available
%> EQuUs v4.8 or later
%%
classdef Custom_Hamiltonians < handle & Messages

properties (Access = protected )
    %> An instance of structure param
    param
    %> Cell array of Hamiltonians of one slab in the leads
    H0 = []; 
    %> Cell array of couplings between the slabs of the leads
    H1 = []; 
    %> Cell array of transverse coupling between the unit cells of the lead
    H1_transverse = [];
    %> Hamiltonian of the scattering region
    Hscatter = []; 
    %> transverse coupling for the scattering region
    Hscatter_transverse = [];
    %> Cell array of couplings between the scattering region and the leads
    Hcoupling = []; 
    %> Cell array of overlap integrals of one slab in the leads
    S0 = []; 
    %> Cell array of overlap integrals between the slabs of the leads
    S1 = []; 
    %> overlap integrals for the transverse coupling between the unit cells of the lead
    S1_transverse = [];
    %> Overlap integrals of the scattering region
    Sscatter = []; 
    %> overlap integrals for the transverse coupling in the scattering region
    Sscatter_transverse = [];
    %> Cell array of overlap integrals between the scattering region and the leads
    Scoupling = []; 
    %> Cell array of coordinates of the leads
    coordinates = []; 
    %> coordinates of the scattering region
    coordinates_scatter = []; 
    %> logical value: true if Hamiltonians are loaded, false otherwise.
    Hamiltonians_loaded = false;
    %> The Fermi energy
    EF
    %> Function Handle for the custom Hamiltonians with identical input values as #LoadHamiltonians, and output values described in the example #Hamiltonians.
    CustomHamiltoniansHandle
    %> list of optional parameters (see #InputParsing for details)
    varargin    
end



methods (Access = public )
    
    
%% Constructor of the class
%> @brief Constructor of the class.
%> @param Opt An instance of the structure #Opt.
%> @param param An instance of structure #param.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
function obj = Custom_Hamiltonians( Opt, param, varargin )
    obj = obj@Messages( Opt );
    obj.param = param;
    obj.varargin = varargin;
    
    obj.Initialize();        
    
end


%% LoadHamiltonians
%> @brief Obtain the Hamiltonians from the external source
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'q' The transverse momentum quantum number.
% @param 'Extmol' The transverse momentum quantum number.
% @param 'GollumInput' The transverse momentum quantum number.
% @param 'WorkingDir' The transverse momentum quantum number.
function LoadHamiltonians(obj, varargin)
    
    p = inputParser;
    
    p.addParameter('q', []);
    p.addParameter('Extmol', 'Extended_Molecule.mat');
    p.addParameter('GollumInput', 'input.mat');
    p.addParameter('WorkingDir', pwd);
    
    p.parse(varargin{:});
    
    
    q = p.Results.q;
    Extmol = p.Results.Extmol;
    GollumInput = p.Results.GollumInput;
    WorkingDir  = p.Results.WorkingDir;    
    
    if isempty(obj.Opt.custom_Hamiltonians)
        % no external source is selected
        warning(['EQuUs:', class(obj), ':LoadHamiltonians'], 'No external source was defined in the input parameters.')
        return
        
    % Siesta Hamiltonians
    elseif strcmpi( obj.Opt.custom_Hamiltonians, 'siesta' )

    
        % Coordinatzes needs to be determined : TODO
        [obj.H0, obj.S0, obj.H1, obj.S1, obj.coordinates, ...
            obj.Hscatter, obj.Sscatter, obj.coordinates_scatter, obj.Hcoupling, obj.Scoupling] = obj.Siesta_Equus_interface(WorkingDir, GollumInput, Extmol);

        %!!!!!!!!!!!!!!!!!! Correct in the Gollum_Equus_interface
        for idx = 1: length(obj.Hcoupling)
            obj.Hcoupling{idx} = obj.Hcoupling{idx}';
            obj.Scoupling{idx} = obj.Scoupling{idx}';            
        end  
    
        obj.H1_transverse = cell(size(obj.H0));
        
        if ~issparse(obj.Hscatter)
           obj.Hscatter = sparse(obj.Hscatter); 
        end
        
        if ~issparse(obj.Sscatter)
           obj.Sscatter = sparse(obj.Sscatter); 
        end     
        
%        obj.S1{1} = obj.S1{1}';
%        obj.S1{2} = obj.S1{2}';
%              
        
        
    %**************************************** develpoment lines *******************************************************
%     obj.H0{2} = obj.H0{1};
%     obj.H1{2} = obj.H1{1};
%     obj.Hscatter = sparse(obj.H0{1});
%     
%     obj.S0{2} = obj.S0{1};
%     obj.S1{2} = obj.S1{1};
%     obj.Sscatter = sparse(obj.S0{1});
    
% obj.H0{1} = obj.H0{1} + 0.5*obj.S0{1};
% obj.H1{1} = obj.H1{1} + 0.5*obj.S1{1};
% obj.H0{2} = obj.H0{2} + speye(size(obj.H0{2}));
%     for idx = 1:length(obj.H0)
%         obj.S0{idx} = [];speye(size(obj.S0{idx}));
%         obj.S1{idx} = [];sparse([],[],[], size(obj.H1{idx},1), size(obj.H1{idx},2));
%     end
%     obj.Sscatter = [];speye(size(obj.Hscatter));
%     for idx =1:length(obj.Hcoupling)
%        obj.Scoupling{idx} = [];sparse([],[],[], size(obj.Hcoupling{idx},1), size(obj.Hcoupling{idx},2));
%     end
    
    
    
    % Custom created Hamiltonians
    elseif strcmpi( obj.Opt.custom_Hamiltonians, 'custom' )
        
        if isempty(obj.CustomHamiltoniansHandle)
            error(['EQuUs:', class(obj), ':LoadHamiltonians'], 'Undefinied function handle for the custom Hamiltonians')
        end
        
        [obj.H0, obj.S0, obj.H1, obj.S1, obj.H1_transverse, obj.coordinates, ...
            obj.Hscatter, obj.Sscatter, obj.Hscatter_transverse, obj.Sscatter_transverse, ...
            obj.coordinates_scatter, obj.Hcoupling, obj.Scoupling] = obj.CustomHamiltoniansHandle();
        
    else
        error(['EQuUs:', class(obj), ':LoadHamiltonians'], 'Undefinied option of the attribute Opt.custom_Hamiltonians')
        
    end
    
    
    % transform Hamiltonians into grand canonical potential
    if ~isempty( obj.EF )        
        for idx = 1:length( obj.S0 )
            if ~isempty( obj.S0{idx} )
                obj.H0{idx} = obj.H0{idx} - obj.EF*obj.S0{idx};
                obj.H1{idx} = obj.H1{idx} - obj.EF*obj.S1{idx};
            else
                obj.H0{idx} = obj.H0{idx} - obj.EF*speye(size(obj.H0{idx}));
            end
        end  
        
        for idx = 1:length( obj.Scoupling )
            if ~isempty( obj.Scoupling{idx} )
                obj.Hcoupling{idx} = obj.Hcoupling{idx} - obj.EF*obj.Scoupling{idx};
            end            
        end
        
        if ~isempty(obj.Sscatter)
            obj.Hscatter = obj.Hscatter - obj.EF*obj.Sscatter;      
        else            
            obj.Hscatter = obj.Hscatter - obj.EF*speye(size(obj.Hscatter));      
        end
    end
    
    
    obj.Hamiltonians_loaded = true;
    
%     for idx = 1:length(obj.H0)
%         obj.S0{idx} = [];speye(size(obj.S0{idx}));
%         obj.S1{idx} = [];sparse([],[],[], size(obj.H1{idx},1), size(obj.H1{idx},2));
%     end
%     obj.Sscatter = [];speye(size(obj.Hscatter));
%     for idx =1:length(obj.Hcoupling)
%        obj.Scoupling{idx} = [];sparse([],[],[], size(obj.Hcoupling{idx},1), size(obj.Hcoupling{idx},2));
%     end    
    
end





%% Reset
%> @brief Resets all elements in the object.
    function Reset( obj )
        
        if strcmpi( class(obj), 'Custom_Hamiltonians' )
            meta_data = metaclass(obj);
            
            for idx = 1:length(meta_data.PropertyList)
                prop_name = meta_data.PropertyList(idx).Name;
                if strcmp(prop_name, 'Opt') || strcmp( prop_name, 'param') || strcmp(prop_name, 'varargin')
                    continue
                end
                obj.Clear( prop_name ); 
            end   
        end 
        
        obj.Initialize();    
    end



%% Write
%> @brief Sets the value of an attribute in the class.
%> @param MemberName The name of the attribute to be set.
%> @param input The value to be set.
    function Write( obj, MemberName, input)
        
        if strcmp(MemberName, 'param')
            obj.param = input;
            obj.Reset();
            return
        elseif strcmp(MemberName, 'Opt')
            obj.Opt = input;
            obj.Reset();           
        else
               try
               obj.(MemberName) = input;
               catch
                    error(['EQuUs:', class(obj), ':Read'], ['No property to write with name: ',  MemberName]);
               end
        end        
    end
%% Read
%> @brief Query for the value of an attribute in the class.
%> @param MemberName The name of the attribute to be set.
%> @return Returns with the value of the attribute.
    function ret = Read(obj, MemberName)
        try
          ret = obj.(MemberName);
          catch
               error(['EQuUs:', class(obj), ':Read'], ['No property to read with name: ',  MemberName]);
          end 
    end
%% Clear
%> @brief Clears the value of an attribute in the class.
%> @param MemberName The name of the attribute to be cleared.
    function Clear(obj, MemberName)
        try
          obj.(MemberName) = [];
          catch
               error(['EQuUs:', class(obj), ':Clear'], ['No property to clear with name: ',  MemberName]);
          end  
        
    end

end % methods public


methods ( Access = private ) 
    
%% Initialize
%> @brief Initializes object properties.
    function Initialize(obj)
        
        obj.Hamiltonians_loaded = false;
        obj.EF = 0;
        
        obj.InputParsing( obj.varargin{:} )        
    end            
    
%% Siesta_Equus_interface
%> @brief Creates the Hamiltonians and overlap integrals from the loaded data from the Siesta package.
%> @param Directory
%> @param input
%> @param extmol
%> @return ...
    function [HL0,SL0,HL1,SL1,coordsLeads,HSC,SSC,coordsSC,HC,SC] = Siesta_Equus_interface( obj, Directory,input,extmol)
        
        %obj.param.scatter

        %SCATTERING REGION
        %----------------------------------------------------------------------------------------
        HSC=0; SSC=0;

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

        HandleForLoad = LoadFromFile( fullfile(Directory, input) );
        atom = HandleForLoad.LoadVariable('atom', 'NoEmpty', 'On');
        leadp = HandleForLoad.LoadVariable('leadp', 'NoEmpty', 'On');
        HandleForLoad.ClearLoadedVariables();

        HandleForLoad = LoadFromFile( fullfile(Directory, extmol) );
        HSM = HandleForLoad.LoadVariable('HSM', 'NoEmpty', 'On');
        iorb = HandleForLoad.LoadVariable('iorb', 'NoEmpty', 'On');
        kpoints_EM = HandleForLoad.LoadVariable('kpoints_EM', 'NoEmpty', 'On');
        nspin = HandleForLoad.LoadVariable('nspin', 'NoEmpty', 'On');
        HandleForLoad.ClearLoadedVariables();
        
        

        numLead=size(leadp,1); 
        iorb = cat(2,atom(iorb(:,1),2:3),iorb);                                      % iorb: Region PL atom-lead atom n l m z
        norb=size(iorb,1);                                                           % Define # of orbitals
        orb_molecule = find(iorb(:,1)==0);                                           % Orbitals in Extended Molecule

        %STRUCTURES FOR COORDINATES
        coordsLeads = cell( 1, length(leadp(:,1)) );
        for i=1:length(leadp(:,1))
            coordsLeads{i} = structures('coordinates');
        end
        coordsSC = structures('coordinates');

        %orb_scissors = intersect(find(iorb(:,1)==0),find(iorb(:,2)==0));                 % Orbitals in Molecule

        %GATHER ORBITALS FOR LEADS
        for lead=1:numLead                                                           % Orbitals in PLs up to TPL
            aux = intersect(find(iorb(:,1)==lead),find(iorb(:,2)==1));
            for PL=2:leadp(lead,1) 
                aux=cat(2,aux,intersect(find(iorb(:,1)==lead),find(iorb(:,2)==PL)));
            end
            switch lead 
                case(1) 
                    orb1=aux; 
                case(2) 
                    orb2=aux; 
                case(3) 
                    orb3=aux; 
                case(4) 
                    orb4=aux; 
            end; 

%    coordsLeads(lead).a=1;
%    if (size(iorb)(2)>7)
%      coordsLeads(lead).x=iorb(aux(:,leadp(lead,2)),8);
%      coordsLeads(lead).y=iorb(aux(:,leadp(lead,2)),9);
%      coordsLeads(lead).z=iorb(aux(:,leadp(lead,2)),10);

%MY WAY OF INTERPRETING THE INPUT FILE
%    coordsLeads(lead).x=iorb(aux(leadp(lead,2)+1:leadp(lead,1)),8);
%    coordsLeads(lead).y=iorb(aux(leadp(lead,2)+1:leadp(lead,1)),9);
%    coordsLeads(lead).z=iorb(aux(leadp(lead,2)+1:leadp(lead,1)),10);
            aux = [];
        end

        %LOAD HAMILTONIANS FROM SEPARATE LEAD FILES
        %-------------------------------------------------------------------------------
        for lead=1:numLead  

    
            if (leadp(lead,3)~=0)      
                HandleForLoad = LoadFromFile( fullfile(Directory, ['Lead_', int2str(lead), '.mat']) );
                HSL          = HandleForLoad.LoadVariable('HSL', 'NoEmpty', 'On');
                kpoints_Lead = HandleForLoad.LoadVariable('kpoints_Lead', 'NoEmpty', 'On');
                HandleForLoad.ClearLoadedVariables();
                switch lead 
                    case(1) 
                        HSL_1=HSL; kpoints_1=kpoints_Lead;
                    case(2) 
                        HSL_2=HSL; kpoints_2=kpoints_Lead;
                    case(3) 
                        HSL_3=HSL; kpoints_3=kpoints_Lead;
                    case(4) 
                        HSL_4=HSL; kpoints_4=kpoints_Lead;
                end;
            end; 
        end



        for spin = 1:nspin  
            for k_EM = 1:size(kpoints_EM,1)
                H = [];
                S = [];
                HL0 = [];
                HL1 = [];
                SL0 = [];
                SL1 = [];
                G = [];
                for i=1:size(HSM,1) 
                    if (find(HSM(i,1)==k_EM))                                         % H & S for given k & spin
                        H(HSM(i,2),HSM(i,3))=HSM(i,4+2*spin)+1.0i*HSM(i,5+2*spin);
                        S(HSM(i,2),HSM(i,3))=HSM(i,4)+1.0i*HSM(i,5);
                    end
                end
%      if ( scissors(1) == 1 )                                             % Add scissor corrections
%         scissor_corrections; 
%      end;                  
                for lead=1:numLead 
                    k_choice = [];
                    kj = [];
                    HSL = [];
                    kpoints_Lead = [];
                    HL01 = [];
                    SL0t = [];
                    HL1t = [];
                    SL1t = [];
                    aux2 = [];
                    
                    if (leadp(lead,3)~=0)
                        switch lead 
                        case(1) 
                            HSL=HSL_1; kpoints_Lead=kpoints_1; 
                        case(2) 
                            HSL=HSL_2; kpoints_Lead=kpoints_2;
                        case(3) 
                            HSL=HSL_3; kpoints_Lead=kpoints_3;
                        case(4) 
                            HSL=HSL_4; kpoints_Lead=kpoints_4;
                        end;
                    
                        for k_Lead=1:size(kpoints_Lead,1)
                            if (abs(kpoints_EM(k_EM,1)+kpoints_Lead(k_Lead,1))<1d-6 && abs(kpoints_EM(k_EM,2)+kpoints_Lead(k_Lead,2))<1d-6) 
                               k_choice=k_Lead; kj=1;
                            elseif (abs(kpoints_EM(k_EM,1)-kpoints_Lead(k_Lead,1))<1d-6 && abs(kpoints_EM(k_EM,2)-kpoints_Lead(k_Lead,2))<1d-6)
                                k_choice=k_Lead; kj=0;
                            end
                        end
                    
                        for i=1:size(HSL,1) 
                            if (find(HSL(i,1)==k_choice))
                                HL0t(HSL(i,2),HSL(i,3))=HSL(i,4+2*spin)           + (-1)^kj*1.0i*HSL(i,5+2*spin);
                                SL0t(HSL(i,2),HSL(i,3))=HSL(i,4)                  + (-1)^kj*1.0i*HSL(i,5);
                                HL1t(HSL(i,2),HSL(i,3))=HSL(i,4+2*(nspin+spin+1)) + (-1)^kj*1.0i*HSL(i,5+2*(nspin+spin+1));
                                SL1t(HSL(i,2),HSL(i,3))=HSL(i,4+2*(nspin+1))      + (-1)^kj*1.0i*HSL(i,5+2*(nspin+1));
                            end;
                            HL0{lead}=HL0t;
                            SL0{lead}=SL0t;
                            HL1{lead}=HL1t;
                            SL1{lead}=SL1t;

 
                        end;
                    
                        aux2 = intersect(find(iorb(:,1)==lead),find(iorb(:,2)==leadp(lead,2)));

%WE HAVE TO APPLY THIS SOMEWHERE (ONLY WHEN USING SEPARATE LEAD CALCULATION)
                    correction(lead)=HL0{lead}(1,1)-H(aux2(1),aux2(1));
                    end; %if
                end; %for numLeads

     %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
                HSC = H(orb_molecule,orb_molecule);
                SSC = S(orb_molecule,orb_molecule);
                p_mol=length(HSC); % Define pointers

                for lead=1:numLead
                    switch lead 
                    case(1) 
                        aux=orb1; 
                    case(2) 
                        aux=orb2; 
                    case(3) 
                        aux=orb3; 
                    case(4) 
                        aux=orb4; 
                    end

     %COUPLING BETWEEN EM AND THE FRIST LEAD UNIT CELL 
                    B=reshape(aux,numel(aux),1); C=aux(:,1);
                    HC{lead}=H(orb_molecule,C);   %Attach coupling to Lead
                    SC{lead}=S(orb_molecule,C);   %Attach coupling to Lead

                    if ( leadp(lead,3) == 0 )
                        HL0{lead}=H(aux(:,leadp(lead,1)),aux(:,leadp(lead,1)));      %K0=<TPL|K|TPL>
                        HL1{lead}=H(aux(:,leadp(lead,1)),aux(:,leadp(lead,1)-1));    %K1=<TPL|K|TPL-1>+coords="in"
                        SL0{lead}=S(aux(:,leadp(lead,1)),aux(:,leadp(lead,1)));        %Needed to normalize states
                        SL1{lead}=S(aux(:,leadp(lead,1)),aux(:,leadp(lead,1)-1));
                    end;   

    
                end;%Lead loop
            end %k loop
        end % spin loop

end %function


%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Optional parameters, see the web documantation.
%> @param 'EF' The Fermi energy.
%> @param 'CustomHamiltoniansHandle' A function handle for the custom Hamiltonians with identical input values as #LoadHamiltonians, and output values described in the example #Hamiltonians.
    function InputParsing( obj, varargin )
    
        p = inputParser;
        p.addParameter('EF', []);
        p.addParameter('CustomHamiltoniansHandle', []);
        
        p.parse(varargin{:});

       
        obj.EF       = p.Results.EF;
        obj.CustomHamiltoniansHandle = p.Results.CustomHamiltoniansHandle;
                
        
    end
    

end % methods private

end