CreateHamiltonians.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - CreateHamiltonians
%    Copyright (C) 2009-2015 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 basic Basic Functionalities
%> @{
%> @file CreateHamiltonians.m
%> @brief A class to create and store Hamiltonian of the scattering region
%> @} 
%> @brief A class to create and store Hamiltonian of the scattering region
%%
classdef CreateHamiltonians < handle & Messages

    properties ( Access = protected )
        %> An instance of structure param
        param
        %> A list of the sites to be kept after decimation.
        kulso_szabfokok
        %> Hscatter - E*Sscatter
        Kscatter
        %> Hscatter_transverse - E*Sscatter_transverse
        Kscatter_transverse        
        %> The matrix of the Hamiltonian.
        Hscatter
        %> The matrix of the Hamiltonian corresponding to the transverse coupling.
        Hscatter_transverse
        %> The matrix of the overlap integrals in the Hamiltonian.
        Sscatter
        %> The matrix of the overlap integrals for the transverse coupling.
        Sscatter_transverse        
        %> The matrix of the Peierls phases.
        fazis_mtx_scatter
        %> The matrix of the Peierls phases for the transverse coupling.
        fazis_mtx_scatter_t
        %> An instance of the structure coordinates.
        coordinates
        %> The length of the scattering region in units of the lattice contant.
        height
        %> The number of sites in the cross section.
        width
        %> The number of sites in one unit cell.
        M
        %> A logical value. True if the Hamiltonian was created, false otherwise.
        HamiltoniansCreated
        %> A logical value. True if the Hamiltonian was decimated, or false otherwise.
        HamiltoniansDecimated
        %> A logical value. True if the overlap integrals were applied, false otherwise.
        OverlapApplied         
        %> A logical value. True if the vector potential was incorporated into the Hamiltonian or false otherwise.
        MagneticFieldApplied
        %> A logical value. True if a gauge transformation was applied on the Hamiltonians, or false otherwise.
        GaugeTransformationApplied
        %> A transverse momentum.
        q
        %> list of optional parameters (see http://www.mathworks.com/help/matlab/ref/varargin.html for details)
        varargin        

    end
    

methods ( Access = public )
%% constructorof 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 = CreateHamiltonians(Opt, param, varargin)
        obj = obj@Messages( Opt );
        obj.Opt = Opt;
        obj.param = param;
        obj.varargin = varargin;

        obj.Initialize();          

end

%% CreateScatterH
%> @brief Creates a Hamiltonian of a rectangle shaped area. The created Hamiltonian and coordinates are stored within the object.
%> @param varargin Cell array of optional parameters:
%> @param 'Scatter_UC' An instance of class #CreateLeadHamiltonians or its subclass
%> @param 'CustomHamiltonian' An instance of class #Custom_Hamiltonians for external Hamiltonian source.
function CreateScatterH( obj, varargin )  
    p = inputParser;
    p.addParameter('Scatter_UC', []);      
    p.addParameter('CustomHamiltonian', []);     
    p.parse(varargin{:});       
    Scatter_UC     = p.Results.Scatter_UC;
    CustomHamiltonian = p.Results.CustomHamiltonian;
        
    shape = obj.param.scatter.shape;
    
    obj.width = shape.width;
    obj.height = shape.height; % number of unit cells!!!!!!!!
    
    if isempty(Scatter_UC)
            
        if ~isempty( obj.Opt.custom_Hamiltonians )
            obj.queryHamiltonians( CustomHamiltonian );   
            return
        elseif strcmpi(obj.Opt.Lattice_Type, 'Square')
            createH = Square_Lead_Hamiltonians();
            [H0,H1,H1_transverse,coordinates_unitcell] = createH.Square_Hamiltonians(obj.param.scatter, obj.width);
            H1_skew_left = [];
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'SSH')
            createH = Square_Lead_Hamiltonians();
            [H0,H1,H1_transverse,coordinates_unitcell] = createH.SSH_Hamiltonians(obj.param.scatter);
            H1_skew_left = [];
            H1_skew_right = [];
        elseif strcmp(obj.Opt.Lattice_Type, 'Lieb')
            createH = Square_Lead_Hamiltonians();
            [H0,H1,H1_transverse,obj.coordinates] = createH.Lieb_Hamiltonians(obj.param.scatter, obj.width);
            H1_skew_left = [];
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'Graphene')
            End_Type = obj.param.scatter.End_Type;
            createH = Hex_Lead_Hamiltonians();
            [H0,H1,H1_transverse,coordinates_unitcell] = createH.Graphene_Hamiltonians(obj.param.scatter, obj.width, End_Type);
            H1_skew_left = [];
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'Graphene_SOC')
            End_Type = obj.param.scatter.End_Type;
            createH = Graphene_SOC_Lead_Hamiltonians();
            [H0, H1, H1_transverse, H1_skew_left, H1_skew_right, coordinates_unitcell] = createH.Graphene_SOC_Hamiltonians(obj.param.scatter, obj.width, End_Type);
        elseif strcmpi(obj.Opt.Lattice_Type, 'Silicene')
            End_Type = obj.param.scatter.End_Type;
            createH = Silicene_Lead_Hamiltonians();
            [H0,H1,H1_transverse,coordinates_unitcell] = createH.Silicene_Hamiltonians(obj.param.scatter, obj.width, End_Type);
            H1_skew_left = [];
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'Graphene_Bilayer')
            End_Type = obj.param.scatter.End_Type;
            createH = Graphene_Bilayer_Lead_Hamiltonians();
            [H0,H1,H1_transverse,coordinates_unitcell] = createH.Graphene_Bilayer_Hamiltonians(obj.param.scatter, obj.width, End_Type);
            H1_skew_left = [];
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'Graphene_Bilayer_2')
            End_Type = obj.param.scatter.End_Type;
            createH = Graphene_Bilayer_Lead_Hamiltonians();
            [H0,H1,H1_transverse,coordinates_unitcell] = createH.Graphene_Bilayer_Hamiltonians2(obj.param.scatter, obj.width, End_Type); 
            H1_skew_left = [];
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'Graphene_Bilayer_3')
            End_Type = obj.param.scatter.End_Type;
            createH = Graphene_Bilayer_Lead_Hamiltonians();
            [H0,H1,H1_transverse,coordinates_unitcell] = createH.Graphene_Bilayer_Hamiltonians3(obj.param.scatter, obj.width, End_Type); 
            H1_skew_left = [];
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'Triangle')
            createH = Triangle_Lead_Hamiltonians();
            [H0, H1, H1_transverse, H1_skew_left, coordinates_unitcell] = createH.Triangle_Hamiltonians(obj.param.scatter, obj.width); 
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'TMDC_Monolayer')
            createH = TMDC_Monolayer_Lead_Hamiltonians();
            [H0, H1, H1_transverse, H1_skew_left, coordinates_unitcell] = createH.TMDC_Monolayer_Hamiltonians(obj.param.scatter, obj.width)
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'TMDC_Monolayer_SOC')
            createH = TMDC_Monolayer_SOC_Lead_Hamiltonians();
            [H0, H1, H1_transverse, H1_skew_left, coordinates_unitcell] = createH.TMDC_Monolayer_SOC_Hamiltonians(obj.param.scatter, obj.width); 
            H1_skew_right = [];
        elseif strcmpi(obj.Opt.Lattice_Type, 'TMDC_Bilayer_SOC')
            createH = TMDC_Bilayer_SOC_Lead_Hamiltonians();
            [H0, H1, H1_transverse, H1_skew_left, coordinates_unitcell] = createH.TMDC_Bilayer_SOC_Hamiltonians(obj.param.scatter, obj.width);
            H1_skew_right = [];
        else
            error('Unrecognized lattice type, or the claculation options are not set to use outher DFT source.')
        end
    else
       H0 = Scatter_UC.Read( 'H0' );
       H1 = Scatter_UC.Read( 'H1' );
       H1_transverse = Scatter_UC.Read( 'H1_transverse' );
       H1_skew_left = Scatter_UC.Read( 'H1_skew_left' );
       H1_skew_right = Scatter_UC.Read( 'H1_skew_right' );
       coordinates_unitcell = Scatter_UC.Read( 'coordinates' ); 
       obj.MagneticFieldApplied  = Scatter_UC.Read( 'MagneticFieldApplied' );
       obj.GaugeTransformationApplied = Scatter_UC.Read( 'GaugeTransformationApplied' );
       obj.width = size(H0, 1);
    end
    
    H1adj = H1';
    obj.M = size(H0,1);
    
    height = obj.height;
    coordinates = structures( 'coordinates' );
    coordinates.a = coordinates_unitcell.a;
    coordinates.b = coordinates_unitcell.b;
    
    % creating coupling between the unit cells
    HCell = repmat({H1}, 1, height-1);
    H1_big = blkdiag(HCell{:});
    H1_big = [ [sparse([],[],[], size(H1_big,1), size(H1,2)), H1_big]; sparse([],[],[], size(H1,1), size(H1_big,2)+size(H1,2))];
    
    % creating coupling between the unit cells in the opposite direction
    HCell = repmat({H1adj}, 1, height-1);
    H1adj_big = blkdiag(HCell{:}); 
    H1adj_big = [sparse([],[],[], size(H1,1), size(H1adj_big,1)+size(H1,2)); [H1adj_big, sparse([],[],[], size(H1adj_big,1), size(H1,2))]]; 
    
    % creating Hamiltonian of the unit cells
    HCell = repmat({H0}, 1, height);
    Hscatter = blkdiag(HCell{:});  
    
    
    % adding coupling between the unit cells       
    Hscatter = Hscatter + H1_big + H1adj_big;
    
    % creating the transverse coupling
    if ~isempty(obj.q)
        HCell = repmat({H1_transverse}, 1, height);
        Hscatter_transverse = blkdiag(HCell{:});
        
        % adding skew_left couplings
        if ~isempty( H1_skew_left )
            HCell = repmat( {H1_skew_left}, 1, height-1 );
            H1_skew_left_big = blkdiag(HCell{:});  
            H1_skew_left_big = [ [sparse([],[],[], size(H1_skew_left_big,1), size(H1_skew_left,2)), H1_skew_left_big]; sparse([],[],[], size(H1_skew_left,1), size(H1_skew_left_big,2)+size(H1_skew_left,2)) ];
            Hscatter_transverse = Hscatter_transverse + H1_skew_left_big;
        end
        
        % adding skew_right couplings
        if ~isempty( H1_skew_right )
            HCell = repmat( {H1_skew_right}, 1, height-1 );
            H1_skew_right_big = blkdiag(HCell{:});  
            H1_skew_right_big = [ sparse([],[],[], size(H1_skew_right,1), size(H1_skew_right_big,2)+size(H1_skew_right,2)); [H1_skew_right_big, sparse([],[],[], size(H1_skew_right_big,1), size(H1_skew_right,2))] ];
            Hscatter_transverse = Hscatter_transverse + H1_skew_right_big;
        end
        
        
                
    else
        Hscatter_transverse = [];
    end
    
    coordfieldnames = fieldnames( coordinates_unitcell );
    for idx = 1:length( coordfieldnames )
        coordfieldname = coordfieldnames{idx};
        if strcmpi( coordfieldname, 'a') || strcmpi( coordfieldname, 'b') || strcmpi( coordfieldname, 'LatticeConstant')
            continue
        elseif strcmpi( coordfieldname, 'x')
            tmp = ones(size(coordinates_unitcell.x))*(coordinates_unitcell.a(1)*(0:height-1));
            tmp = reshape(tmp, numel(tmp), 1);
            coordinates.x = repmat(coordinates_unitcell.x, height, 1) + tmp;
        elseif strcmpi( coordfieldname, 'y')
            tmp = ones(size(coordinates_unitcell.y))*(coordinates_unitcell.a(2)*(0:height-1));
            tmp = reshape(tmp, numel(tmp), 1);    
            coordinates.y = repmat(coordinates_unitcell.y, height, 1) + tmp;
        elseif ~isempty(coordinates_unitcell.(coordfieldname))
            coordinates.(coordfieldname) = repmat(coordinates_unitcell.(coordfieldname), height, 1);
        end
    end
    
    
    obj.coordinates = coordinates;
    obj.Hscatter = Hscatter;
    obj.Hscatter_transverse = Hscatter_transverse;
    
    obj.HamiltoniansCreated   = true;
    obj.MagneticFieldApplied  = false;
    obj.GaugeTransformationApplied = false;
    obj.OverlapApplied = false;
    obj.HamiltoniansDecimated = false;
    
    
end


%% ApplyOverlapMatrices
%> @brief Applies the overlap matrices to the Hamiltonians: K = H-ES 
%> @param E The energy value.
    function ApplyOverlapMatrices(obj, E)
        
        if obj.OverlapApplied
            obj.display(['EQuUs:', class(obj), ':ApplyOverlapMatrices: Overlap matrices were already applied.']);
            return;
        end
        
        if ~isempty( obj.Sscatter )
            obj.Kscatter = obj.Hscatter - E*obj.Sscatter;     
            obj.OverlapApplied = true;
        else
            obj.Kscatter = obj.Hscatter - eye(size(obj.H0))*E;
        end
                
        
        if ~isempty( obj.Sscatter_transverse )
            obj.Kscatter_transverse = obj.Hscatter_transverse - E*obj.Sscatter_transverse;
        elseif ~isempty( obj.H1_transverse )
            obj.Kscatter_transverse = obj.Hscatter_transverse;
        end        
                        
    end

%% shiftFermiEnergy 
%> @brief Shifts the on-site energies in the Scattering region by a given energy.
%> @param Energy The energy value to be used.
    function shiftFermiEnergy( obj, Energy )  
        if isempty( obj.Sscatter )
            energyshift = speye( size( obj.Hscatter ) )*Energy;
        else
            energyshift = obj.Sscatter*Energy;
        end        
        
        obj.Hscatter = obj.Hscatter - energyshift;
        if ~isempty( obj.param.scatter) && ~isempty( obj.param.scatter.epsilon )
            obj.param.scatter.epsilon = obj.param.scatter.epsilon + Energy;
        end
    end

%% projectHamiltonian2Spin
%> @brief Projects the Hamiltonian to a spin states of $s=\pm1$
%> @param s The quantum index of the spin (\pm1)
    function projectHamiltonian2Spin( obj, s )
        if s == 1
            obj.coordinates.x = obj.coordinates.x(obj.coordinates.spinup);
            obj.coordinates.y = obj.coordinates.y(obj.coordinates.spinup);
            obj.Hscatter = obj.Hscatter( obj.coordinates.spinup, obj.coordinates.spinup);
            if ~isempty(obj.q)
                obj.Hscatter_transverse = obj.Hscatter_transverse( obj.coordinates.spinup, obj.coordinates.spinup);
            end
        elseif s == -1
            obj.coordinates.x = obj.coordinates.x(~obj.coordinates.spinup);
            obj.coordinates.y =obj. coordinates.y(~obj.coordinates.spinup);
            obj.Hscatter = obj.Hscatter( ~obj.coordinates.spinup, ~obj.coordinates.spinup);    
            if ~isempty(obj.q)
                obj.Hscatter_transverse = obj.Hscatter_transverse( ~obj.coordinates.spinup, ~obj.coordinates.spinup);
            end
        else
            display(['EQuUs:', class(obj), ':projectHamiltonian2Spin: spin should be +/-1'], 1)
        end
        
        obj.coordinates.spinup = [];
        obj.M = obj.M/2;
        
        
    end
    
%% Reset
%> @brief Resets all elements in the class.
    function Reset(obj)
        
        if strcmpi( class(obj), 'CreateHamiltonians' )
            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

%% CreateClone
%> @brief Creates a clone of the present class.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'empty' Set true to create an empty clone, or false (default) to clone all atributes.
%> @return Returns with the cloned object.
    function ret = CreateClone( obj, varargin )

        p = inputParser;
        p.addParameter('empty', false ); 
        
        p.parse(varargin{:});       
        empty    = p.Results.empty;
        
        ret = CreateHamiltonians(obj.Opt, obj.param, obj.varargin{:});
                                        
        if empty
            return
        end
                  
        meta_data = metaclass(obj);
        
        for idx = 1:length(meta_data.PropertyList)
            prop_name = meta_data.PropertyList(idx).Name;
               ret.Write( prop_name, obj.(prop_name));  
        end   
                      
    end   
    
    
%% RemoveSites
%> @brief Removes specific sites from the model of the scattering region. For example see function #ScatterPotQD.
%> @param indexes Logical array. Sites with true values will be removed from the system.
    function RemoveSites( obj, indexes )
        
        % update the array of external degrees of freedom
        if ~isempty( obj.kulso_szabfokok )
            all_sites = 1:size(obj.Hscatter,1);
            sites2remove = all_sites( indexes );
            
            for idx = length(sites2remove):-1:1
               indexes_tmp = find( obj.kulso_szabfokok == sites2remove(idx), 1 );
               if ~isempty(indexes_tmp)
                   error(['EQuUs:', class(obj), ':RemoveSites'], 'Sites cannot be removed from the Hamiltonian.');
               end
               
               indexes_tmp = obj.kulso_szabfokok > sites2remove(idx);
               obj.kulso_szabfokok(indexes_tmp) = obj.kulso_szabfokok(indexes_tmp) - 1;
            end
        end
        
        % removing sites from the scattering Hamiltonian
        obj.Hscatter = obj.Hscatter(~indexes, ~indexes);
        
        % removing sites from the matrix of the transverse coupling
        if ~isempty(obj.Hscatter_transverse)
            obj.Hscatter_transverse = obj.Hscatter_transverse(~indexes, ~indexes);
        end
        
        % removing sites from the matrix of the overlap integrals
        if ~isempty(obj.Sscatter )
            obj.Sscatter = obj.Sscatter(~indexes, ~indexes);
        end
        
        % removing sites from the matrix of the overlap integrals of the transverse coupling
        if ~isempty(obj.Sscatter_transverse )
            obj.Sscatter_transverse = obj.Sscatter_transverse(~indexes, ~indexes);
        end
        
        % removing sites from the matrix of the Peierls phases
        if ~isempty(obj.fazis_mtx_scatter )
            obj.fazis_mtx_scatter = obj.fazis_mtx_scatter(~indexes, ~indexes);
        end
        
        % removing sites from the matrix of the Peierls phases of the transverse coupling
        if ~isempty(obj.fazis_mtx_scatter_t )
            obj.fazis_mtx_scatter_t = obj.fazis_mtx_scatter_t(~indexes, ~indexes);
        end
        
        % removing sites from structure coordinates
        obj.coordinates = obj.coordinates.RemoveSites( indexes );        

        
    end       

%% saveScatterHamiltonian
%> @brief Saves the Hamiltonian and other data of the scattering region
%> @param filename The string containing the path to the file. (In octave use absolute paths only)
    function saveScatter( filename )
        if ~exist('filename', 'var')
            filename = 'Hscatter.mat';
        end
        save(filename,'Hscatter', 'Hscatter_transverse', 'coordinates', 'param', 'fazis_mtx_scatte', 'fazis_mtx_scatte_t', 'width', 'height', 'M')
    end

%% LoadHamiltonian
%> @brief Loads the Hamiltonian and other data of the scattering region from a file.
%> @param filename The string containing the path to the file. (In octave use absolute paths only)
    function loadScatter( obj, filename )
    if ~exist('filename', 'var')
        filename = 'Hscatter.mat';
    end    
    disp('EQuUs::CreateHamiltonian::loadScatter: Loading scatter Hamiltonian from file')
    HandleForLoad = LoadFromFile( filename );
    obj.Hscatter = HandleForLoad.LoadVariable('Hscatter', 'NoEmpty', 'On');
    obj.Hscatter_transverse = HandleForLoad.LoadVariable('Hscatter_transverse', 'NoEmpty', 'On');
    obj.Sscatter = HandleForLoad.LoadVariable('Sscatter', 'NoEmpty', 'On');
    obj.Sscatter_transverse = HandleForLoad.LoadVariable('Sscatter_transverse', 'NoEmpty', 'On');    
    obj.coordinates = HandleForLoad.LoadVariable('coordinates', 'NoEmpty', 'On');
    obj.fazis_mtx_scatter = HandleForLoad.LoadVariable('fazis_mtx_scatte', 'NoEmpty', 'Off');
    obj.fazis_mtx_scatter_t = HandleForLoad.LoadVariable('fazis_mtx_scatte_t', 'NoEmpty', 'Off');
    obj.param = HandleForLoad.LoadVariable('param', 'NoEmpty', 'On');
    obj.M = HandleForLoad.LoadVariable('M', 'NoEmpty', 'On');
    obj.width = HandleForLoad.LoadVariable('width', 'NoEmpty', 'Off');
    obj.height = HandleForLoad.LoadVariable('height', 'NoEmpty', 'On');
     obj.q = HandleForLoad.LoadVariable('q', 'NoEmpty', 'On');
    HandleForLoad.ClearLoadedVariable();
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, 'params')
            obj.param.scatter = input;
            return            
        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

%% MemberClear
%> @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 attributes.
    function Initialize(obj)
        
        obj.coordinates          = structures('coordinates');
        obj.HamiltoniansCreated   = false; 
        obj.HamiltoniansDecimated = false;
        obj.MagneticFieldApplied  = false;
        obj.GaugeTransformationApplied = false;
        obj.OverlapApplied = false;   
        
        obj.InputParsing( obj.varargin{:} );        
    end        
    
%% queryHamiltonians
%> @brief Obtains Hamiltonians from a class instance Custom_Hamiltonians
%> @param cCustom_Hamiltonians An instance of class #Custom_Hamiltonians
    function queryHamiltonians( obj, cCustom_Hamiltonians ) 
        if isempty( cCustom_Hamiltonians )
            cCustom_Hamiltonians = Custom_Hamiltonians( obj.Opt, obj.param ); 
        end
        
        if ~cCustom_Hamiltonians.Read( 'Hamiltonians_loaded' )
            cCustom_Hamiltonians.LoadHamiltonians();
        end
            
        obj.coordinates         = cCustom_Hamiltonians.Read( 'coordinates_scatter' );
        obj.Hscatter            = cCustom_Hamiltonians.Read( 'Hscatter' );
        obj.Hscatter_transverse = cCustom_Hamiltonians.Read( 'Hscatter_transverse' );
        obj.Sscatter            = cCustom_Hamiltonians.Read( 'Sscatter' );
        obj.Sscatter_transverse = cCustom_Hamiltonians.Read( 'Sscatter_transverse' );        
        obj.HamiltoniansCreated   = true;
        obj.MagneticFieldApplied  = false;
        obj.GaugeTransformationApplied = false;  
        obj.HamiltoniansDecimated = false;
        
        % determine kulso_szabfokok
        Hcoupling = cCustom_Hamiltonians.Read( 'Hcoupling' );
        
        if iscell( Hcoupling )
            for idx = 1:length( Hcoupling )
               coupled_sites = CoupledSites( Hcoupling{idx} );
               obj.kulso_szabfokok = [obj.kulso_szabfokok, coupled_sites];
            end
            obj.kulso_szabfokok = unique( obj.kulso_szabfokok );
        else
            obj.kulso_szabfokok = CoupledSites( Hcoupling );
        end
        
        
        % tranforming into the BdG model if necessary
        if ~isempty( obj.Opt.BdG ) && obj.Opt.BdG
            obj.Transforms2BdG()
        end
        
        % nested functions
        
        % determine sites coupled to the leads
        function ret = CoupledSites( Hcoupling )
           [rows, cols] =  find(Hcoupling);
           ret = unique( cols );
           ret = reshape(ret, 1, numel(ret));
        end
        
        % end nested functions
    end
    
%% Transforms2BdG
%> @brief Transforms the Hamiltonians into a Bogoliubov de Gennes model.
    function Transforms2BdG( obj )
        if isempty( obj.param.scatter.pair_potential )
            error('CreateHamiltonian:Transform2BdG', 'Pair potential need to be set for the scattering region.');
        end
        
        if ~isempty(obj.coordinates.BdG_u)
            obj.display('Hamiltonians already transformed to BdG model.')
            return
        end
        
        
        
        % transforming the Hamiltonians
        pair_potential = obj.param.scatter.pair_potential;        
        site_num = size( obj.Hscatter, 1 );        
        obj.coordinates.BdG_u = [true( site_num, 1); false( site_num, 1) ];
        
        if isempty( obj.Sscatter )
            Sscatter = speye(size(obj.Hscatter));
        else
            Sscatter = obj.Sscatter;
        end
        
        Hscatter_hole = -conj(obj.Hscatter);
        % apply time reversal symmetry on spins according to Eq (3) of PRB 86 134522 (2012)
        if ~isempty(obj.coordinates.spinup)
            Hscatter_hole = [Hscatter_hole(~obj.coordinates.spinup, ~obj.coordinates.spinup), -Hscatter_hole(~obj.coordinates.spinup, obj.coordinates.spinup);
                             -Hscatter_hole(obj.coordinates.spinup, ~obj.coordinates.spinup), Hscatter_hole(obj.coordinates.spinup, obj.coordinates.spinup) ];
        end
        obj.Hscatter = [ obj.Hscatter, Sscatter*pair_potential;
                         Sscatter*conj(pair_potential), Hscatter_hole ];
        
        % transforming the transverse coupling
        if ~isempty(obj.Hscatter_transverse)    
            Hscatter_transverse_hole = -conj(obj.Hscatter_transverse);
            
            % apply time reversal symmetry on spins according to Eq (3) of PRB 86 134522 (2012)
            if ~isempty(obj.coordinates.spinup)
                Hscatter_transverse_hole = [Hscatter_transverse_hole(~obj.coordinates.spinup, ~obj.coordinates.spinup), -Hscatter_transverse_hole(~obj.coordinates.spinup, obj.coordinates.spinup);
                             -Hscatter_transverse_hole(obj.coordinates.spinup, ~obj.coordinates.spinup), Hscatter_transverse_hole(obj.coordinates.spinup, obj.coordinates.spinup) ];
            end
        
            obj.Hscatter_transverse = [ obj.Hscatter_transverse, sparse([],[],[], size(obj.Hscatter_transverse,1), size(obj.Hscatter_transverse,2));
                         sparse([],[],[], size(obj.Hscatter_transverse,1), size(obj.Hscatter_transverse,2)), Hscatter_transverse_hole ];
        end
        
        % transforming the overlap integrals
        if ~isempty( obj.Sscatter )
            obj.Sscatter = [ obj.Sscatter, sparse([],[],[], site_num, site_num );
                         sparse([],[],[], site_num, site_num ), obj.Sscatter ];
        end
                     
                     
        fnames = fieldnames( obj.coordinates );
        for idx = 1:length(fnames)
            fname = fnames{idx};
            if strcmp(fname, 'a') || strcmp(fname, 'b')
                continue
            end                    
            obj.coordinates.(fname) = [ obj.coordinates.(fname); obj.coordinates.(fname) ];
        end
                
        obj.coordinates.BdG_u = [ true(size(obj.Hscatter,1),1), false(size(obj.Hscatter,1),1)];
                
        length_kulso_szabfokok = length(obj.kulso_szabfokok);
        obj.kulso_szabfokok = [obj.kulso_szabfokok(1:length_kulso_szabfokok), obj.kulso_szabfokok(1:length_kulso_szabfokok)+size(obj.Hscatter,1) ];          
        
    end
    
%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Optional parameters, see the web documantation.
%> @param 'q' The transverse momentum quantum number.
    function InputParsing( obj, varargin)
        
        p = inputParser;
        p.addParameter('q', []);
        
        p.parse(varargin{:});

        obj.q = p.Results.q;

    end


end %methods protected

end % Global end