2 % Copyright (C) 2016 Peter Rakyta, Ph.D.
4 % This program is free software: you can redistribute it and/or modify
5 % it under the terms of the GNU General Public License as published by
6 % the Free Software Foundation, either version 3 of the License, or
7 % (at your option) any later version.
9 % This program is distributed in the hope that it will be useful,
10 % but WITHOUT ANY WARRANTY; without even the implied warranty of
11 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 % GNU General Public License
for more details.
14 % You should have received a copy of the GNU General Public License
15 % along with
this program. If not, see http:
17 %> @addtogroup basic Basic Functionalities
20 %> @brief Class to regulerize the H1 problem of the leads by SVD decompozition.
22 %> @brief Class to regulerize the H1 problem of the leads by SVD decompozition.
23 %> The notations and the structure of the Hamiltonian is defined accroding to the following image:
24 %> @image html Lead_Hamiltonian.jpg
25 %> @image latex Lead_Hamiltonian.jpg
27 %> EQuUs v4.8 or later
32 properties ( Access =
protected )
33 %>
true if the
Hamiltonians were SVD transformed,
false otherwise
35 %> U matrix from the SVD decompozition, see Eq (41) of PRB 78, 035407
37 %> S matrix from the SVD decompozition, see Eq (41) of PRB 78, 035407
39 %> V matrix from the SVD decompozition, see Eq (41) of PRB 78, 035407
41 %> SVD tolerance to identify singular values
43 %> Somethimes it is needed to perform another SVD cycle to regularize the H1 matrix
45 %> Effective number of
sites after the elimination of the singular values
55 methods ( Access = public )
56 %% constructorof the class
57 %> @brief Constructor of the class.
58 %> @
param Opt An instance of the structure
#Opt. 61 %> @
return An instance of the
class 69 %> @brief Calculates the SVD decomposition of the matrix H1.
70 function SVDdecompozition( obj )
72 [obj.U, obj.S, obj.V] = svd(full(obj.K1));
74 [obj.U, obj.S, obj.V] = svd(obj.K1);
77 obj.U = sparse(obj.U);
78 obj.V = sparse(obj.V);
84 %> @brief Decides whether SVD regularization is needed or not.
85 %> @
return Returns with
true if SVD regularization is needed,
false otherwise
86 function ret = is_SVD_needed(obj)
87 if 1/condest(obj.H1)<obj.tolerance
95 %% Calc_Effective_Hamiltonians
96 %> @brief Calculates the effective
Hamiltonians according to Eq (48) of of PRB 78, 035407
97 %> @
param E The energy value
98 function Calc_Effective_Hamiltonians( obj, E )
100 obj.ApplyOverlapMatrices(E);
102 if ~isempty(obj.kulso_szabfokok) && length(obj.kulso_szabfokok) ~= size(obj.K0,1)
103 obj.Decimate_Hamiltonians();
105 elseif 1/condest(obj.K1)<obj.tolerance && isempty(obj.V)
109 % in case no SVD regularization or simple decimation is needed
110 obj.Neff = size(obj.K0,1);
118 %% Get_Effective_Hamiltonians
119 %> @brief Gets the effective
Hamiltonians K0_eff, K1_eff, K1adj_eff according to Eq (48) of of PRB 78, 035407
120 %> @return [1] The effective Hamiltonian of one unit cell
121 %> @return [2] The effective coupling between unit cells K1_eff
122 %> @return [3] The adjungate of the effective coupling between unit cells K1adj_eff
123 function [K0_eff, K1_eff, K1adj_eff] = Get_Effective_Hamiltonians(obj)
124 if isempty(obj.next_SVD_cycle)
127 K1adj_eff = obj.K1adj;
129 [K0_eff, K1_eff, K1adj_eff] = obj.next_SVD_cycle.Get_Effective_Hamiltonians();
131 error('EQuUs:
SVDregularizationLead:Get_Effective_Hamiltonians', 'Unrecognized type of attribute next_SVD_cycle');
136 %% Get_Effective_Overlaps
137 %> @brief Gets the effective
Hamiltonians S0_eff, S1_eff, S1adj_eff according to Eq (48) of of PRB 78, 035407
138 %> @return [1] The effective overlap matrix of one unit cell S0_eff
139 %> @return [2] The effective overlap matrix between unit cells S1_eff
140 %> @return [3] The adjungate of the effective overlap matrix between unit cells S1adj_eff
141 function [S0_eff, S1_eff, S1adj_eff] = Get_Effective_Overlaps(obj)
142 if isempty(obj.next_SVD_cycle)
145 S1adj_eff = obj.S1adj;
147 [S0_eff, S1_eff, S1adj_eff] = obj.next_SVD_cycle.Get_Effective_Overlaps();
149 error('EQuUs:
SVDregularizationLead:Get_Effective_Overlaps', 'Unrecognized type of attribute next_SVD_cycle');
154 %% Get_Effective_Coordinates
155 %> @brief Gets the coordinates of the
sites of the effective
Hamiltonians. (Has sense if the singular
sites were given directly)
156 %> @return [1] A class
Coordinates containing the coordinates of the
sites of the effective Hamiltonian
157 function coordinates = Get_Effective_Coordinates(obj)
158 if isempty(obj.next_SVD_cycle)
159 coordinates = obj.coordinates;
161 coordinates = obj.next_SVD_cycle.Get_Effective_Coordinates();
163 error('EQuUs:
SVDregularizationLead:Get_Effective_Coordinates', 'Unrecognized type of attribute next_SVD_cycle');
170 %> @brief Regularize the
Hamiltonians of the lead by SVD regularization
173 obj.SVDdecompozition();
177 if length(obj.q) < 2 && ~isempty(obj.q) && ~obj.HamiltoniansDecimated
178 K0 = K0 + obj.K1_transverse*diag(exp(1i*obj.q)) + obj.K1_transverse'*diag(exp(-1i*obj.q));
181 % Eqs (48) in PRB 78, 035407
182 % here we use different transformation for the leads with different orientation
183 if obj.Lead_Orientation == 1
185 elseif obj.Lead_Orientation == -1
190 K1adj = U'*obj.K1adj*U;
192 obj.Neff = size(obj.K0, 1);
194 % determine singular values
196 regular_sites_slab = transpose(find(S > obj.tolerance*max(S)));
198 K = [K0, K1; K1adj, K0];
201 regular_sites = [regular_sites_slab, size(K0,1)+regular_sites_slab];
204 CreateH.Write('Hscatter', K);
205 CreateH.Write('kulso_szabfokok', regular_sites);
206 CreateH.Write('HamiltoniansCreated', true);
207 CreateH.Write('HamiltoniansDecimated', false);
210 Decimation_Procedure.DecimationFunc(0, CreateH, 'Hscatter', 'kulso_szabfokok');
212 K = CreateH.Read('Hscatter');
213 CreateH.Clear('Hscatter');
215 Neff_new = length(regular_sites_slab);
217 if obj.Lead_Orientation == 1
218 K0 = K(Neff_new+1:end, Neff_new+1:end);
219 elseif obj.Lead_Orientation == -1
220 K0 = K(1:Neff_new, 1:Neff_new);
222 K1 = K(1:Neff_new, Neff_new+1:end);
223 K1adj = K(Neff_new+1:end, 1:Neff_new);
228 obj.next_SVD_cycle.Write(
'K1', K1);
229 obj.next_SVD_cycle.Write(
'K1adj', K1adj);
230 obj.next_SVD_cycle.Write(
'OverlapApplied',
true);
231 obj.next_SVD_cycle.Write(
'HamiltoniansDecimated',
true);
232 obj.next_SVD_cycle.Write(
'Neff', Neff_new);
235 obj.is_SVD_transformed =
true;
236 obj.HamiltoniansDecimated =
true;
239 condition_num = 1/condest(K1);
241 condition_num = rcond(K1);
244 if condition_num<obj.tolerance
245 obj.next_SVD_cycle.Calc_Effective_Hamiltonians( 0 );
249 %% Decimate_Hamiltonians
251 function Decimate_Hamiltonians( obj )
253 K0loc = obj.K0 + obj.K1_transverse*diag(exp(1i*obj.q)) + obj.K1_transverse'*diag(exp(-1i*obj.q));
257 K = [K0loc,obj.K1;obj.K1adj,K0loc];
259 regular_sites_slab = obj.kulso_szabfokok;
260 regular_sites = [regular_sites_slab, size(obj.K0,1)+regular_sites_slab];
264 CreateH.Write('Hscatter', K);
265 CreateH.Write('kulso_szabfokok', regular_sites);
266 CreateH.Write('HamiltoniansCreated', true);
267 CreateH.Write('HamiltoniansDecimated', false);
270 Decimation_Procedure.DecimationFunc(0, CreateH, 'Hscatter', 'kulso_szabfokok');
272 K = CreateH.Read('Hscatter');
273 CreateH.Clear('Hscatter');
274 coordinates = CreateH.Read('coordinates');
276 Neff_new = length(regular_sites_slab);
278 K0 = K(Neff_new+1:end, Neff_new+1:end);
279 K1 = K(1:Neff_new, Neff_new+1:end);
280 K1adj = K(Neff_new+1:end, 1:Neff_new);
284 obj.next_SVD_cycle.Write(
'K1', K1);
285 obj.next_SVD_cycle.Write(
'K1adj', K1adj);
286 obj.next_SVD_cycle.Write(
'OverlapApplied',
true);
287 obj.next_SVD_cycle.Write(
'HamiltoniansDecimated',
true);
288 obj.next_SVD_cycle.Write(
'Neff', Neff_new);
291 obj.HamiltoniansDecimated =
true;
294 % determine the truncated coordinates
295 indexes2remove =
true( size( obj.coordinates.x ) );
296 indexes2remove( regular_sites_slab ) =
false;
297 coordinates = obj.coordinates.RemoveSites( indexes2remove );
298 obj.next_SVD_cycle.Write(
'coordinates', coordinates);
302 % Decimation_Procedure_Lead =
Decimation(obj.Opt,
'ForcedDecimation', obj.Opt.Decimation);
303 % kulso_szabfokok_tmp = obj.kulso_szabfokok;
304 % obj.kulso_szabfokok = [obj.kulso_szabfokok, obj.kulso_szabfokok+size(obj.H0,1)];
305 % Decimation_Procedure_Lead.DecimationFunc(0, obj,
'Hamiltonian2Dec',
'kulso_szabfokok');
306 % obj.kulso_szabfokok = kulso_szabfokok_tmp;
307 % obj.Neff = size(obj.K0,1);
312 %> @brief Transforms the effective
Hamiltonians by a unitary transformation
313 %> @
param Umtx The matrix of the unitary transformation.
314 function Unitary_Transform(obj, Umtx)
316 if isempty(obj.next_SVD_cycle)
318 obj.K0 = Umtx*obj.K0*Umtx';
322 obj.K1 = Umtx*obj.K1*Umtx';
325 if ~isempty(obj.K1adj)
326 obj.K1adj = Umtx*obj.K1adj*Umtx';
330 obj.next_SVD_cycle.Unitary_Transform(Umtx);
332 error('EQuUs:
SVDregularizationLead:Unitary_Transform', 'Unrecognized type of attribute next_SVD_cycle');
338 %> @brief Gets the effective number of
sites after the elimination of the singular values.
339 %> @return Returns with the effective number of
sites 341 if isempty(obj.next_SVD_cycle)
344 Neff = obj.next_SVD_cycle.Get_Neff();
354 %> @brief Gets the total transformation U related to the SVD transformation
355 %> @return Returns with the total transformation U
357 if isempty(obj.next_SVD_cycle)
360 Vtot_tmp = obj.next_SVD_cycle.Get_V();
361 size_V = size(obj.V);
362 size_Vtmp = size(Vtot_tmp);
363 Vtot_tmp = [Vtot_tmp, zeros(size_Vtmp(1), size_V(2)-size_Vtmp(2));
364 zeros(size_V(1)-size_Vtmp(1), size_Vtmp(2)), eye(size_V(1)-size_Vtmp(1))];
365 Vtot = obj.V*Vtot_tmp;
373 %> @brief Creates a clone of the present
object.
374 %> @return Returns with the cloned
object.
375 %> @
param varargin Cell array of optional parameters:.
376 function ret = CreateClone( obj, varargin )
379 p.addParameter('empty', false );
381 p.parse(varargin{:});
382 empty = p.Results.empty;
385 'Hanyadik_Lead', obj.Hanyadik_Lead,...
386 'Lead_Orientation', obj.Lead_Orientation, ...
393 meta_data = metaclass(obj);
395 for idx = 1:length(meta_data.PropertyList)
396 prop_name = meta_data.PropertyList(idx).Name;
397 if strcmpi( prop_name, 'next_SVD_cycle' )
398 if ~isempty( obj.next_SVD_cycle )
399 ret.next_SVD_cycle = obj.next_SVD_cycle.CreateClone();
402 ret.Write( prop_name, obj.(prop_name));
410 %> @brief Resets all elements in the
object.
414 meta_data = metaclass(obj);
416 for idx = 1:length(meta_data.PropertyList)
417 prop_name = meta_data.PropertyList(idx).Name;
418 if strcmp(prop_name, '
Opt') || strcmp( prop_name, '
param') || strcmp(prop_name, 'varargin')
421 obj.Clear( prop_name );
432 %> @brief Sets the value of an attribute in the interface.
433 %> @
param MemberName The name of the attribute to be set.
434 %> @
param input The value to be set.
435 function Write(obj, MemberName, input)
438 if strcmpi(MemberName, '
param') || strcmp(MemberName, 'params')
443 obj.(MemberName) = input;
445 error(['EQuUs:', class(obj), ':Read'], ['No property to write with name: ', MemberName]);
451 %> @brief Query for the value of an attribute in the interface.
452 %> @
param MemberName The name of the attribute to be set.
453 %> @return Returns with the value of the attribute.
454 function ret = Read(obj, MemberName)
456 if strcmpi(MemberName, 'Hamiltonian2Dec')
461 ret = obj.(MemberName);
463 error(['EQuUs:', class(obj), ':Read'], ['No property to read with name: ', MemberName]);
469 %> @brief Clears the value of an attribute in the interface.
470 %> @
param MemberName The name of the attribute to be cleared.
474 obj.(MemberName) = [];
476 error(['EQuUs:', class(obj), ':Clear'], ['No property to clear with name: ', MemberName]);
486 methods (Access=protected)
489 %% Extend_Wavefnc ---- DEVELOPMENT
490 %> @brief Extend a reduced wave
function to the original basis before the SVD regularization (Eq (45) in PRB 78 035407
491 %> @
param wavefnc_reduced The reduced wavefunction of the effective system
492 %> @
param expk e^(i*k)
493 %> @return A matrix conatining the extended wave functions.
494 function wavefnc_extended = Extend_Wavefnc(obj, wavefnc_reduced, expk)
497 wavefnc_reduced = obj.next_SVD_cycle.Extend_Wavefnc(wavefnc_reduced, expk);
498 elseif ~isempty(obj.next_SVD_cycle)
499 error('EQuUs:
SVDregularizationLead:Extend_Wavefnc', 'Unrecognized type of attribute next_SVD_cycle');
502 Fn = -obj.invD*(obj.K1adju/expk + obj.C);
503 wavefnc_extended = [wavefnc_reduced; Fn*wavefnc_reduced];
504 wavefnc_extended = obj.U*wavefnc_extended;
511 %> @brief Initializes class properties.
514 obj.tolerance = 1e-12;
515 obj.is_SVD_transformed = false;
521 end % methods protected
Property next_SVD_cycle
Somethimes it is needed to perform another SVD cycle to regularize the H1 matrix.
Structure Opt contains the basic computational parameters used in EQuUs.
Class to create and store Hamiltonian of the translational invariant leads.
function Transport(Energy, B)
Calculates the conductance at a given energy value.
function Hamiltonians(varargin)
Function to create the custom Hamiltonians for the 1D chain.
A class providing function handle to reduce the number of sites in the Hamiltonian via decimation pro...
Class to regulerize the H1 problem of the leads by SVD decompozition.
function CreateLeadHamiltonians(Opt, param, varargin)
Constructor of the class.
Structure param contains data structures describing the physical parameters of the scattering center ...
Structure sites contains data to identify the individual sites in a matrix.
Structure containing the coordinates and other quantum number identifiers of the sites in the Hamilto...
A class to create and store Hamiltonian of the scattering region.