2 % Copyright (C) 2009-2015 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 utilities Utilities
20 %> @brief A
class to perform transport calculations on a graphene
antidot (i.e., a hollow in a ribbon). Obsolete
class,
for real calculations a creation of a
new class is recommended.
22 %> @brief A
class to perform transport calculations on a graphene
antidot (i.e., a hollow in a ribbon). Obsolete class, for real calculations a creation of a new class is recommended.
26 properties ( Access =
protected )
27 %> A lattice vector in the hexagonal lattice.
29 %> A lattice vector in the hexagonal lattice.
31 %> Directory name to export the plotted wave functions.
33 %> The name of the subdirectory to export the plotted wave functions.
37 properties ( Access =
public )
38 %> An instance of structure
hole 42 %> The strength of the magnetic field in Tesla
48 %> flux in the
hole in units of \phi_0
54 %> The charge of the electron
60 methods ( Access =
public )
63 %> @brief Constructor of the
class.
65 %> @
param interpq An array of transverse momentum points at which fhandle is calculated via interpolation.
66 %> @
param varargin Cell array of optional parameters. For details see #InputParsing.
67 %> @
return An instance of the
class 68 function obj =
antidot( varargin )
70 obj.
a1 = [sqrt(3); -1]*sqrt(3)/2;
71 obj.a2 = [sqrt(3); 1]*sqrt(3)/2;
77 obj.B = 10e-10; %in Tesla
78 obj.antidot_edge_points = [];
79 obj.flux_of_hole = [];
85 obj.waveFncDirnameFull = [];
86 obj.waveFncSubDirname = [];
88 obj.InputParsing( varargin{:} );
89 obj.generate_geometry();
91 obj.display(
'Creating antidot object')
103 %> @brief Calculates the conductance of an
antidot connected to the leads in the
"contact/scattering center/contact" arrangement.
104 %> @
param Energy The energy value.
105 %> @
return C The calculated conductivity
106 %> @
return ny Array of the open channel in the leads.
107 %> @
return DeltaC Error of the unitarity.
108 function [C,ny,DeltaC] =
Transport( obj, Energy )
110 obj.display(
'Calculating transport in create_Hole_in_ribbon')
111 obj.setEnergy( Energy );
112 obj.create_scatter_GreensFunction();
114 Dysonfunc = @()obj.CustomDysonFunc(
'gfininv', obj.Ginv,
'SelfEnergy',
false,
'constant_channels',
false);
116 obj.FL_handles.DysonEq(
'CustomDyson', Dysonfunc );
118 save(
'Error_antidot_Transport.mat');
119 for idx = 1:length(errCause.stack)
120 obj.display(errCause.stack(idx), 1)
124 [S,ny] = obj.FL_handles.SmatrixCalc();
126 norma = norm(S*S'-eye(sum(ny)));
128 obj.display( ['error of the unitarity of S-matrix: ',num2str(norma)] )
132 obj.display( ['openchannels do not match: ',num2str(ny)] )
135 conductance = obj.FL_handles.Conduktance();
136 conductance = abs([conductance(1,:), conductance(1,:)]);
137 DeltaC = std(conductance);
139 C = mean(conductance);
141 obj.display( ['Energy = ', num2str(Energy), ' conductance = ', num2str(C)])
147 %> @brief Calculates the energies and wave functions of the bound states.
148 %> @
param Energy The energy value.
149 %> @
param closefigure If true, the created figure is closed at the end of the calculations.
150 %> @
param varargin Cell array of optional parameters:
151 %> @
param 'toPlot' Set true (default) for plotting the wave function, or false otherwise.
152 %> @
param 'Einhomegac' Set true if the Energy is given in units of $$\hbar\omega_c$$, or false (default) otherwise.
153 %> @
param 'filterHole' Set true (default) for separating
antidot bound states from the edge states, or false otherwise.
154 %> @
param 'db' The number of the calculated energy eigenvalues. Default is db=30.
155 %> @
param 'infinitemass' Set true for using the infinite mass boundary condition, or false (default) otherwise.
156 %> @
param 'dotCalc' Set true for a dot, or false (default) otherwise.
157 %> @
param 'delta' Every delta-th point of the wave function becomes plotted. Default is delta=4.
158 %> @
param 'smoothedge' Set true for using a smooth edge in the calculations, or false (default) otherwise.
159 %> @return ret A structure containing the calculated results.
160 %TODO create structure corresponding to the return value
161 function ret = CalcWavefnc( obj, Energy, closefigure, varargin )
162 obj.display( 'Calculating the wave function' );
163 if ~exist('closefigure', 'var')
168 p.addParameter('toPlot', 1, @isscalar);
169 p.addParameter('Einhomegac', 0, @isscalar);
170 p.addParameter('filterHole', 1, @isscalar);
171 p.addParameter('db', 30, @isscalar);
172 p.addParameter('infinitemass', 0, @isscalar);
173 p.addParameter('dotCalc', 0, @isscalar);
174 p.addParameter('delta', 4, @isscalar);
175 p.addParameter('smoothedge', 0, @isscalar);
177 p.parse(varargin{:});
178 toPlot = p.Results.toPlot;
179 Einhomegac = p.Results.Einhomegac;
180 filterHole = p.Results.filterHole;
182 infinitemass = p.Results.infinitemass;
183 dotCalc = p.Results.dotCalc;
184 delta = p.Results.delta; %used in wavefunction plotting
185 smoothedge = p.Results.smoothedge;
188 rCC = obj.rCC*1e-10; %Angstrom
189 phi0 = obj.h/obj.qe; % magnetic flux quantum
190 eta_B = 2*pi/phi0*(rCC)^2*obj.B;
191 homega = 2.97*sqrt(9/2*eta_B);
192 Energy = Energy*homega;
197 if dotCalc && ~infinitemass
200 [Hscatter, wavefunc_indexes] = obj.create_hole_Hamiltonian(
'infinitemass', infinitemass,
'dotCalc', dotCalc,
'smoothedge', smoothedge );
201 xcoords = obj.coordinates.x
'; 202 ycoords = obj.coordinates.y';
205 height = obj.height*2;
207 obj.display(
' Calculating eigenvalues and wave functions')
208 [eigvecs,eigvals] = eigs(Hscatter,db,Energy);
210 eigvals = diag( eigvals );
211 [eigvals, IX] = sort(eigvals);
212 eigvecs = eigvecs(:,IX);
213 ret =
struct(
'eigenval',[],
'Wavefnc', [],
'xcoords', [],
'ycoords', []);
219 obj.display( eigvals)
220 ret.eigenval = eigvals;
222 % reshaping the eigenvectors
223 Wavefnc = cell( size(eigvals ) );
224 for kdx = 1:length(eigvals)
226 eigvec = eigvecs(:,kdx);
227 wavefnc = zeros(size(xcoords));
229 % plotting wavefunction calculated by the faster way
230 sorok = round( (wavefunc_indexes - mod(wavefunc_indexes,M))/M ) + 1;
231 indexes_enrows = logical( mod(wavefunc_indexes,M) == 0 );
232 sorok( indexes_enrows ) = sorok( indexes_enrows ) -1;
235 indexes_in_row = logical( sorok == idx );
236 indexes_in_row = wavefunc_indexes(indexes_in_row) - M*(idx-1);
237 wavefnc(idx, indexes_in_row ) = eigvec( indexek_tmp + (1:length(indexes_in_row) ) );
238 indexek_tmp = indexek_tmp + length(indexes_in_row);
241 Wavefnc{kdx} = wavefnc;
244 ret.eigenval = eigvals;
245 ret.Wavefnc = Wavefnc;
246 ret.xcoords = xcoords;
247 ret.ycoords = ycoords;
255 mkdir( obj.waveFncDirnameFull );
256 mkdir( obj.waveFncDirnameFull, obj.waveFncSubDirname);
257 obj.display(
' Plotting wavefunctions' )
258 for kdx = 1:length(eigvals)
261 figure1 = figure(
'XVisual',...
262 '0x21 (TrueColor, depth 24, RGB mask 0xff0000 0xff00 0x00ff)',...
263 'Renderer',
'painters',...
264 'Colormap',[1 1 1;0.948412716388702 0.948412716388702 0.948412716388702;0.896825432777405 0.896825432777405 0.896825432777405;0.845238089561462 0.845238089561462 0.845238089561462;0.793650805950165 0.793650805950165 0.793650805950165;0.742063522338867 0.742063522338867 0.742063522338867;0.690476179122925 0.690476179122925 0.690476179122925;0.638888895511627 0.638888895511627 0.638888895511627;0.58730161190033 0.58730161190033 0.58730161190033;0.581477105617523 0.581477105617523 0.581477105617523;0.575652658939362 0.575652658939362 0.575652658939362;0.569828152656555 0.569828152656555 0.569828152656555;0.564003705978394 0.564003705978394 0.564003705978394;0.558179199695587 0.558179199695587 0.558179199695587;0.552354753017426 0.552354753017426 0.552354753017426;0.546530246734619 0.546530246734619 0.546530246734619;0.540705800056458 0.540705800056458 0.540705800056458;0.534881293773651 0.534881293773651 0.534881293773651;0.52905684709549 0.52905684709549 0.52905684709549;0.523232340812683 0.523232340812683 0.523232340812683;0.517407834529877 0.517407834529877 0.517407834529877;0.511583387851715 0.511583387851715 0.511583387851715;0.505758881568909 0.505758881568909 0.505758881568909;0.499934434890747 0.499934434890747 0.499934434890747;0.494109958410263 0.494109958410263 0.494109958410263;0.488285452127457 0.488285452127457 0.488285452127457;0.482460975646973 0.482460975646973 0.482460975646973;0.476636499166489 0.476636499166489 0.476636499166489;0.470812022686005 0.470812022686005 0.470812022686005;0.464987546205521 0.464987546205521 0.464987546205521;0.459163069725037 0.459163069725037 0.459163069725037;0.445249050855637 0.445249050855637 0.445249050855637;0.431335002183914 0.431335002183914 0.431335002183914;0.417420983314514 0.417420983314514 0.417420983314514;0.403506934642792 0.403506934642792 0.403506934642792;0.389592915773392 0.389592915773392 0.389592915773392;0.375678867101669 0.375678867101669 0.375678867101669;0.361764848232269 0.361764848232269 0.361764848232269;0.347850799560547 0.347850799560547 0.347850799560547;0.333936780691147 0.333936780691147 0.333936780691147;0.320022732019424 0.320022732019424 0.320022732019424;0.306108713150024 0.306108713150024 0.306108713150024;0.292194694280624 0.292194694280624 0.292194694280624;0.278280645608902 0.278280645608902 0.278280645608902;0.264366626739502 0.264366626739502 0.264366626739502;0.25045257806778 0.25045257806778 0.25045257806778;0.236538544297218 0.236538544297218 0.236538544297218;0.222624525427818 0.222624525427818 0.222624525427818;0.208710491657257 0.208710491657257 0.208710491657257;0.194796457886696 0.194796457886696 0.194796457886696;0.180882424116135 0.180882424116135 0.180882424116135;0.166968390345573 0.166968390345573 0.166968390345573;0.153054356575012 0.153054356575012 0.153054356575012;0.139140322804451 0.139140322804451 0.139140322804451;0.12522628903389 0.12522628903389 0.12522628903389;0.111312262713909 0.111312262713909 0.111312262713909;0.0973982289433479 0.0973982289433479 0.0973982289433479;0.0834841951727867 0.0834841951727867 0.0834841951727867;0.0695701614022255 0.0695701614022255 0.0695701614022255;0.0556561313569546 0.0556561313569546 0.0556561313569546;0.0417420975863934 0.0417420975863934 0.0417420975863934;0.0278280656784773 0.0278280656784773 0.0278280656784773;0.0139140328392386 0.0139140328392386 0.0139140328392386;0 0 0]);
269 axes1 = axes(
'Parent',figure1,
'CLim',[0 0.028]);
272 contourf( xcoords(1:delta:end,1:delta:end), ycoords(1:delta:end,1:delta:end), abs(wavefnc(1:delta:end,1:delta:end)), 50,
'LineStyle',
'none',
'Parent', axes1 );
275 text(30,30, [
'E = ', num2str(eigvals(kdx)*1000),
' meV']);
276 set(gcf,
'Renderer',
'painters')
278 filename = ['wavefunc_width',num2str(obj.width),'_height',num2str(obj.height),'_radius',num2str(obj.
hole.radius),'_B',num2str(obj.B),'_E',num2str(eigvals(kdx),'%6.5f'),'.eps'];
280 filename = ['infmass_',filename];
283 filename = ['smoothedge_',filename];
286 filename = ['dot_', filename];
288 print('-depsc2', [obj.waveFncDirnameFull, filesep, obj.waveFncSubDirname,'/',filename]);
300 %----------------------------------------------
303 %----------------------------------
304 %this filters out eigenstates related to bound states to the Hole
305 function filterHoleStates()
306 % disp( ' Filtering out
hole states' )
307 hole_indexes = true(1,length(eigvals));
310 for iidx = 1:length(hole_indexes)
311 if norm(eigvecs(1:M-20,iidx)) > tol
312 hole_indexes(iidx) = false;
316 eigvals = eigvals( hole_indexes );
317 eigvecs = eigvecs(:, hole_indexes );
320 % end nested functions
325 %% create_hole_Hamiltonian
326 %> @brief Creates the Hamiltonian for the
antidot/dot.
327 %> @
param varargin Cell array of optional parameters:
328 %> @
param 'infinitemass' Set true for using the infinite mass boundary condition, or false (default) otherwise.
329 %> @
param 'dotCalc' Set true for a dot, or false (default) otherwise.
330 %> @
param 'smoothedge' Set true for using a smooth edge in the calculations, or false (default) otherwise.
331 %> @return Hscatter The matrix representation of the created Hamiltonian.
332 %> @return wavefunction_indexes The list of the
sites included in the
antidot/dot system.
333 function [Hscatter, wavefunction_indexes] = create_hole_Hamiltonian( obj, varargin )
335 p.addParameter('infinitemass', 0, @isscalar);
336 p.addParameter('dotCalc', 0, @isscalar);
337 p.addParameter('smoothedge', 0, @isscalar);
338 p.parse(varargin{:});
340 infinitemass = p.Results.infinitemass;
341 dotCalc = p.Results.dotCalc;
342 smoothedge = p.Results.smoothedge;
347 obj.display(
' Creating the Hamiltonian of a finite ribbon')
348 obj.CreateScatter( );
349 Hscatter = obj.CreateH.Read(
'Hscatter' );
351 obj.display(
' Creating hole in the ribbon ')
355 Hscatter_tmp = Hscatter;
358 Hscatter_tmp = Hscatter;
361 Hscatter_tmp = Hscatter;
366 if isempty(obj.antidot_edge_points)
370 obj.display( ' cutting out
hole points ')
373 for idx = size(obj.antidot_edge_points,1):-1:1
374 indexek_tmp = (obj.antidot_edge_points(idx)-1)*(2*M) + (obj.antidot_edge_points(idx,2) : obj.antidot_edge_points(idx,3));
375 Hscatter( indexek_tmp,: ) = 0;
376 Hscatter( :, indexek_tmp ) = 0;
377 indexes2clear = [indexes2clear, indexek_tmp];
380 indexes2clear_tmp = getHoleSitesFromTheOddRows( indexes2clear );
381 indexes2clear = sort( [ indexes2clear, indexes2clear_tmp]);
385 indexes_tmp = 1:size(Hscatter_tmp,1);
386 indexes_tmp( indexes2clear) = [];
389 wavefunction_indexes = indexes2clear;
390 Hscatter = Hscatter_tmp;
391 Hscatter( indexes_tmp,: ) = [];
392 Hscatter( :, indexes_tmp ) = [];
395 wavefunction_indexes = indexes_tmp;
396 Hscatter( indexes2clear,: ) = [];
397 Hscatter( :, indexes2clear ) = [];
399 wavefunction_indexes = 1:size(Hscatter,1);
404 indexes_tmp = 1:size(Hscatter_tmp,1);
405 indexes_tmp( indexes2clear ) = [];
406 indexes2clear = indexes_tmp;
408 ApplyInfiniteMass( indexes2clear )
409 wavefunction_indexes = 1:size(Hscatter_tmp,1);
412 %----------------------------------------------
415 %-----------------------------------------
416 function ApplyInfiniteMass( indexek)
417 Hscatter = Hscatter_tmp;
419 sorok = round((indexek - mod(indexek,M))/M)+1;
420 indexes_endrows = logical( mod(indexek,M) == 0 );
421 sorok( indexes_endrows ) = sorok( indexes_endrows ) -1;
423 ps_sorok = logical( mod(sorok,2) == 0 );
424 pt_sorok = ~ps_sorok;
426 ps_sorok = indexek( ps_sorok );
427 pt_sorok = indexek( pt_sorok );
429 idx_A = [ ps_sorok( logical( mod(ps_sorok,2) == 0)), pt_sorok( logical( mod(pt_sorok,2) == 1)) ];
430 idx_B = [ ps_sorok( logical( mod(ps_sorok,2) == 1)), pt_sorok( logical( mod(pt_sorok,2) == 0)) ];
432 Hscatter = Hscatter + sparse( idx_A, idx_A, V, size(Hscatter,1), size(Hscatter,2)) + sparse( idx_B, idx_B, -V, size(Hscatter,1), size(Hscatter,2));
438 % ----------------------------------------
439 function indexes2clear = getHoleSitesFromTheOddRows( indexek_tmp )
443 for iidx = 1:length( indexek_tmp )
444 related_indexes = find(Hscatter_tmp(indexek_tmp(iidx),:));
445 sorok = round((related_indexes - mod(related_indexes,M))/M)+1;
446 indexes_endrows = logical( mod(related_indexes,M) == 0 );
447 sorok( indexes_endrows ) = sorok( indexes_endrows ) -1;
449 ps_sorok = logical( mod(sorok,2) == 0);
450 related_indexes = related_indexes( ps_sorok );
452 for iiidx = 1:length(related_indexes)
453 idx_tmp = find(Hscatter(related_indexes(iiidx),:));
454 idx_tmp2 = find( idx_tmp ~= related_indexes(iiidx) );
455 if length( idx_tmp2 ) <= 1;
456 if isempty( find(indexes2clear == related_indexes(iiidx),1) )
457 indexes2clear = [indexes2clear, related_indexes(iiidx)];
467 % end nested functions
473 %> @brief Creates function handles of the vector potentials and apply the magnetic filed in the ribbon
Hamiltonians.
476 obj.setHandlesForMagneticField('scatter', hLandauy, 'lead', hLandauy );
482 methods ( Access = protected )
485 %> @brief creates
handle for vector potential
488 rCC = obj.rCC*1e-10; %Angstrom
491 eta_B = 2*pi/phi0*(rCC)^2*B;
492 Aconst = 0.1; % constant vektor potential in the systems: stabilizes the numerical computation
494 hLandauy = @(x,y)([zeros(size(x,1),1),eta_B*x]);
502 %> @brief creates geometry data of the
hole in the ribbon
503 function generate_geometry( obj )
506 obj.display(' Generating geometry in
antidot ')
509 % creating geometry with the
510 obj.convert_site_indexes( )
513 end % methods private
515 methods ( Access = public )
518 %> @brief Determines the
sites that should be cut off from the ribbon in order to create the
hole.
519 %> @returns antidot_edge_points A matrix with culomns z (slab index), zpoints_min, zpoints_max (lower and upper bounds of the
sites in the slab z).
520 function antidot_edge_points = getHolePoints( obj )
522 if isempty(obj.Scatter_UC)
523 obj.CreateRibbon('justHamiltonians', true);
526 M = obj.Scatter_UC.Read( 'M' );
528 deltax_ps = 2; % 2n-1 es 2n kozotti tavolsag rCC egysegekben
529 deltax_pt = 1; % 2n es 2n+1 kozotti tavolsag rCC egysegekben
530 deltax = 3; % 2n-1 es 2n+1 kozotti tavolsag rCC egysegekben
533 hole.center = [obj.
hole.center(1)*(deltax_ps+deltax_pt)/2; obj.
hole.center(2)*deltay];
536 hole.radius = obj.
hole.radius*(norm(rA2))/2;
538 zmin = round((
hole.center(2) -
hole.radius)/deltay) - 1;
539 zmax = round((
hole.center(2) +
hole.radius)/deltay) + 1;
541 antidot_edge_points = zeros( zmax-zmin, 3);
547 dist = distance_from_center(ztmp, idx );
548 if ( dist<
hole.radius && ~minfound )
553 if ( dist>
hole.radius && minfound )
555 antidot_edge_points(hole_edge_idx,:) = [ztmp, idx_min, idx_max];
556 hole_edge_idx = hole_edge_idx + 1;
563 antidot_edge_points(hole_edge_idx:end,:) = [];
564 obj.antidot_edge_points = antidot_edge_points;
566 %----------------------------------------------
569 %----------------------------------------------
570 function ret = distance_from_center(z, idx )
572 maradek = mod(idx,2);
575 xcoord = (idx-1)/2*deltax;
577 xcoord = (idx-2)/2*deltax + deltax_ps;
581 ret = sqrt( (xcoord-
hole.center(1))^2 + (ycoord -
hole.center(2))^2 );
585 % end nested functions
591 %% create_scatter_GreensFunction
592 %> @brief Calculates the surface Greens function of the
antidot. Sets the attributes
#gfin and #gfininv to the calculated results. 593 function create_scatter_GreensFunction( obj )
595 obj.Scatter_UC.ApplyOverlapMatrices(obj.E);
596 K0 = obj.Scatter_UC.Read(
'K0');
597 K1 = obj.Scatter_UC.Read(
'K1');
598 K1adj = obj.Scatter_UC.Read(
'K1adj');
602 obj.display([
'EQuUs:Utils:',
class(obj),
':CalcFiniteGreensFunction: Eigenvalues of the scattering region'])
603 obj.Scatter_UC.TrukkosSajatertekek(obj.E);
605 obj.display([
'EQuUs:Utils:',
class(obj),
':CalcFiniteGreensFunction: Group velocity for the scattering region'])
606 obj.Scatter_UC.Group_Velocity();
610 if isempty(obj.antidot_edge_points)
614 %> the first two and last two slabs are added manualy at the end
619 obj.G = GreensFuncAtRibbonEnds();
620 antidot_edge_points = obj.antidot_edge_points;
621 scatterer_points = getScattererPoints();
623 if ~isempty( scatterer_points )
624 ghole = GreensFuncAtHoleEdge();
625 [ghopp1, ghopp2] = GreensFuncAtCouplings();
626 obj.G = [obj.G,ghopp1; ghopp2,ghole];
629 if isnan(rcond(obj.G) ) || abs( rcond(obj.G) ) < 1e-15
630 obj.Ginv= obj.inv_SVD( obj.G );
632 obj.Ginv= inv(obj.G);%gfin\eye(size(gfin));%inv(gfin);
637 %applying scattering potentials
638 for ldx = 1:length(obj.
scatterers.siteindexes)
642 obj.Ginv= obj.Ginv( [1+M:3*M,obj.
scatterers.siteindexes+4*M], [1+M:3*M,obj.
scatterers.siteindexes+4*M] );
643 kulso_szabfokok = 1:2*M;
644 obj.Ginv= obj.DecimationFunction( kulso_szabfokok, gfininv);
646 obj.Ginv= obj.Ginv( 1+M:3*M, 1+M:3*M );
649 % adding to the scattering region the first and last slabs
650 Neff = obj.Scatter_UC.Get_Neff();
651 non_singular_sites = obj.Scatter_UC.Read( 'kulso_szabfokok' );
652 if isempty( non_singular_sites )
653 non_singular_sites = 1:Neff;
656 % getting the Hamiltonian for the edge slabs
657 non_singular_sites_edges = [non_singular_sites, size(K0,1)+1:2*size(K0,1)];
658 ginv_edges = -[K0, K1; K1adj, K0];
659 ginv_edges = obj.DecimationFunction( non_singular_sites_edges, ginv_edges );
661 %Ginv = [first_slab, H1, 0;
663 % 0 , H1', last_slab];
665 obj.Ginv = [ginv_edges(1:Neff,1:Neff), [ginv_edges(1:Neff, Neff+non_singular_sites), zeros(Neff, Neff+size(K0,2))]; ...
666 [ginv_edges(Neff+non_singular_sites, 1:Neff); zeros(Neff, Neff)], obj.Ginv, [zeros(Neff, size(ginv_edges,2)-Neff); ginv_edges(1:Neff, Neff+1:end)]; ...
667 [zeros(size(K0,2),2*Neff), ginv_edges(Neff+1:end, 1:Neff)], ginv_edges(Neff+1:end, Neff+1:end)];
669 [rows, cols] = find( K1 );
670 rows = unique(rows); % cols identical to non_singular_sites
672 %non_singular_sites_Ginv = [1:length(non_singular_sites), size(obj.Ginv,1)-size(K0,1)+1:size(obj.Ginv,1)];
673 non_singular_sites_Ginv = [1:length(non_singular_sites), size(obj.Ginv,1)-size(K0,1)+reshape(rows, 1, length(rows))];
674 obj.Ginv = obj.DecimationFunction( non_singular_sites_Ginv, obj.Ginv );
676 % apply gauge transformation if given
677 coordinates_scatter = obj.getCoordinates();
680 % gauge transformation on Green's function
681 if ~isempty(obj.Ginv) && ~isempty(obj.PeierlsTransform_Scatter) && isempty(obj.q)
682 % gauge transformation on the inverse Green's function
683 obj.Ginv= obj.PeierlsTransform_Scatter.gaugeTransformation( obj.Ginv, coordinates_scatter, obj.
gauge_field );
686 err = MException('
Ribbon:CalcFiniteGreensFunction', 'Unable to perform gauge transformation');
687 err = addCause(err, errCause);
688 save('Error_Ribbon_CalcFiniteGreensFunction.mat')
693 %----------------------------------------------
697 %--------------------------------------------------------------------------
698 % get
sites of
hole and scatterer that should be obtained in the Green's function
699 function scatterer_points = getScattererPoints()
701 scatterer_points = cell(size(antidot_edge_points,1),3);
703 if ~isempty(antidot_edge_points)
704 point_num_down = antidot_edge_points(1,3)-antidot_edge_points(1,2) + 1;
705 point_num_up = antidot_edge_points(end,3)-antidot_edge_points(end,2) + 1;
708 for idx = 1:size( antidot_edge_points,1 )
709 ztmp = antidot_edge_points( idx );
712 zpoints = antidot_edge_points(idx,2):antidot_edge_points(idx,3);
713 indexek = 1:point_num_down;
714 elseif idx == size(antidot_edge_points,1)
715 zpoints = antidot_edge_points(idx,2):antidot_edge_points(idx,3);
716 indexek = max(indexek)+1:max(indexek)+point_num_up;
718 if antidot_edge_points(idx,2) + 6 < antidot_edge_points(idx,3)
719 zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,2)+2,antidot_edge_points(idx,2)+3,antidot_edge_points(idx,3)-3,antidot_edge_points(idx,3)-2,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
720 indexek = max(indexek)+1:max(indexek)+8;
721 elseif antidot_edge_points(idx,2) + 5 < antidot_edge_points(idx,3)
722 zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,2)+2,antidot_edge_points(idx,2)+3,antidot_edge_points(idx,3)-2,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
723 indexek = max(indexek)+1:max(indexek)+7;
724 elseif antidot_edge_points(idx,2) + 4 < antidot_edge_points(idx,3)
725 zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,2)+2,antidot_edge_points(idx,3)-2,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
726 indexek = max(indexek)+1:max(indexek)+6;
727 elseif antidot_edge_points(idx,2) + 3 < antidot_edge_points(idx,3)
728 zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,2)+2,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
729 indexek = max(indexek)+1:max(indexek)+5;
730 elseif antidot_edge_points(idx,2) + 2 < antidot_edge_points(idx,3)
731 zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,3)-1,antidot_edge_points(idx,3)];
732 indexek = max(indexek)+1:max(indexek)+4;
733 elseif antidot_edge_points(idx,2) + 1 < antidot_edge_points(idx,3)
734 zpoints = [antidot_edge_points(idx,2),antidot_edge_points(idx,2)+1,antidot_edge_points(idx,3)];
735 indexek = max(indexek)+1:max(indexek)+3;
737 zpoints = antidot_edge_points(idx,2:3);
738 indexek = max(indexek)+1:max(indexek)+2;
743 scatterer_points{idx,1} = antidot_edge_points(idx,1);
744 scatterer_points{idx,2} = zpoints;
745 scatterer_points{idx,3} = ones(size(zpoints)); %1
for hole edge points, 0 otherwise
750 z = antidot_edge_points(:,1);
752 for idx = 1:length( obj.scatterers.z )
753 index = find( z == round(obj.
scatterers.z(idx)/2),1 );
755 is_hole_edge = 0; %1 for
hole edge points, 0 otherwise
756 scatterer_points = [scatterer_points; {round(obj.scatterers.z(idx)/2), obj.scatterers.zpoints(idx), is_hole_edge} ];
757 z = [z; round(obj.scatterers.z(idx)/2)];
759 scatterer_points{index,2} = [scatterer_points{index,2}, obj.scatterers.zpoints(idx)];
760 is_hole_edge = 0; %1
for hole edge points, 0 otherwise
761 scatterer_points{index,3} = [scatterer_points{index,3}, is_hole_edge];
765 getScattererSiteIndexes();
769 %-----------------------------------------------------------
770 function getScattererSiteIndexes()
775 for iidx = 1:size(scatterer_points,1)
777 scatterer_idx = find( scatterer_points{iidx,3}==0 );
778 if ~isempty( scatterer_idx )
779 obj.scatterers.
siteindexes(kkdx+1:kkdx+length(scatterer_idx) ) = index_sum + scatterer_idx;
780 kkdx = kkdx + length(scatterer_idx);
782 index_sum = index_sum + length(scatterer_points{iidx,3});
790 %-------------------------------------------------
791 % calculate the Greens
function for the
hole edge points
792 function ghole = GreensFuncAtHoleEdge( )
795 for idx = 1:size(scatterer_points,1)
796 point_num = point_num + length(scatterer_points{idx,2});
799 ghole = zeros(point_num,point_num);
801 for idx = 1:size( scatterer_points,1 )
802 z1tmp = scatterer_points{idx,1};
803 z1points = scatterer_points{idx,2};
804 indexek1 = max(indexek1)+1:max(indexek1)+length(z1points);
807 for jdx = 1:size( scatterer_points,1 )
808 z2tmp = scatterer_points{jdx,1};
809 z2points = scatterer_points{jdx,2};
810 indexek2 = max(indexek2)+1:max(indexek2)+length(z2points);
813 gtmp = obj.Scatter_UC.InfGreenFunction(z1tmp, z2tmp,
'z1points', z1points,
'z2points', z2points);
815 ghole(indexek1, indexek2) = gtmp;
819 ghole(max(indexek1)+1:end,:) = [];
820 ghole(:,max(indexek2)+1:end) = [];
825 %-------------------------------------------------
826 % calculate the Greens
function between the ribbon and points and
hole edge
827 function [ghopp1, ghopp2] = GreensFuncAtCouplings()
830 point_num = size( ghole,1 );
831 M = obj.Scatter_UC.Read( 'M' );
832 ghopp1 = zeros(4*M,point_num);
833 ghopp2 = zeros(point_num,4*M);
836 for idx = 1:size( scatterer_points,1 )
837 ztmp = scatterer_points{idx,1};
838 zpoints = scatterer_points{idx,2};
839 indexek = max(indexek)+1:max(indexek)+length(zpoints);
841 gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, ztmp,
'z2points', zpoints);
842 ghopp1(1:M, indexek) = gtmp;
844 gtmp = obj.Scatter_UC.InfGreenFunction(z1, ztmp,
'z2points', zpoints);
845 ghopp1(1+M:2*M, indexek) = gtmp;
847 gtmp = obj.Scatter_UC.InfGreenFunction(z2, ztmp,
'z2points', zpoints);
848 ghopp1(1+2*M:3*M, indexek) = gtmp;
850 gtmp = obj.Scatter_UC.InfGreenFunction(z2+1, ztmp,
'z2points', zpoints);
851 ghopp1(1+3*M:4*M, indexek) = gtmp;
856 gtmp = obj.Scatter_UC.InfGreenFunction( ztmp, z1-1,
'z1points', zpoints);
857 ghopp2(indexek,1:M) = gtmp;
859 gtmp = obj.Scatter_UC.InfGreenFunction( ztmp, z1,
'z1points', zpoints);
860 ghopp2(indexek,1+M:2*M) = gtmp;
862 gtmp = obj.Scatter_UC.InfGreenFunction( ztmp, z2,
'z1points', zpoints);
863 ghopp2(indexek,1+2*M:3*M) = gtmp;
865 gtmp = obj.Scatter_UC.InfGreenFunction( ztmp, z2+1,
'z1points', zpoints);
866 ghopp2(indexek,1+3*M:4*M) = gtmp;
872 %-------------------------------------------------
873 function gfin = GreensFuncAtRibbonEnds()
875 M = obj.Scatter_UC.Read('M');
876 gfin = zeros(4*M,4*M);
878 gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, z1-1);
879 gfin(1:M,1:M) = gtmp;
880 gfin(1+M:2*M, 1+M:2*M) = gtmp;
881 gfin(1+2*M:3*M, 1+2*M:3*M) = gtmp;
882 gfin(1+3*M:4*M, 1+3*M:4*M) = gtmp;
884 gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, z1);
885 gfin(1:M,1+M:2*M) = gtmp;
886 gfin(1+2*M:3*M,1+3*M:4*M) = gtmp;
888 gtmp = obj.Scatter_UC.InfGreenFunction(z1, z1-1);
889 gfin(1+M:2*M,1:M) = gtmp;
890 gfin(1+3*M:4*M, 1+2*M:3*M) = gtmp;
892 gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, z2);
893 gfin(1:M,1+2*M:3*M) = gtmp;
895 gtmp = obj.Scatter_UC.InfGreenFunction(z2, z1-1);
896 gfin(1+2*M:3*M,1:M) = gtmp;
898 gtmp = obj.Scatter_UC.InfGreenFunction(z1-1, z2+1);
899 gfin(1:M,1+3*M:4*M) = gtmp;
901 gtmp = obj.Scatter_UC.InfGreenFunction(z2+1, z1-1);
902 gfin(1+3*M:4*M,1:M) = gtmp;
904 gtmp = obj.Scatter_UC.InfGreenFunction(z1, z2);
905 gfin(1+M:2*M,1+2*M:3*M) = gtmp;
907 gtmp = obj.Scatter_UC.InfGreenFunction(z2, z1);
908 gfin(1+2*M:3*M,1+M:2*M) = gtmp;
910 gtmp = obj.Scatter_UC.InfGreenFunction(z1, z2+1);
911 gfin(1+M:2*M,1+3*M:4*M) = gtmp;
913 gtmp = obj.Scatter_UC.InfGreenFunction(z2+1, z1);
914 gfin(1+3*M:4*M,1+M:2*M) = gtmp;
918 % end nested functions
927 methods ( Access = private )
929 %% convert_site_indexes
930 %> @brief Converets site indexes into coordinates
931 function convert_site_indexes( obj )
933 xcoord = zeros( obj.width, obj.height*2 );
934 ycoord = zeros( obj.width, obj.height*2 );
941 radius = obj.
hole.radius*(norm(rA2))/2;
943 obj.flux_of_hole = pi*radius^2*(obj.rCC*1e-10)^2*obj.B/(obj.h/obj.qe);
948 for idx = 1:obj.width
949 for jdx = 1:obj.height*2
951 if mod(idx+jdx,2) == 0 % A atom
952 r =
double(idx-1)/2*rA1 +
double(jdx-1)/2*rA2;
954 r =
double(idx)/2*rA1 +
double(jdx-1)/2*rA2 + rB;
957 xcoord(idx,jdx) = r(1);
958 ycoord(idx,jdx) = r(2);
964 %center = obj.
hole.center(1)/2*rA1 + obj.
hole.center(2)/2*rA2;
965 %obj.coordinates.outhole = logical( (xcoord - center(1)).^2 + (ycoord - center(2)).^2 >= radius^2 );
967 obj.coordinates.x = xcoord*obj.rCC;
968 obj.coordinates.y = ycoord*obj.rCC;
975 end % methods private
977 methods ( Access = protected )
980 %> @brief Parses the optional parameters for the class constructor.
981 %> @
param varargin Cell array of optional parameters:
982 %> @
param 'width' Integer. The number of the atomic
sites in the cross section of the ribbon.
983 %> @
param 'height' Integer. The height of the ribbon in units of the lattice vector.
984 %> @
param 'filenameIn' Input filename for the xml input structure.
985 %> @
param 'filenameOut' Output filename for the xml input structure.
986 %> @
param 'WorkingDir' The absolute path to the working directoy.
987 %> @
param 'E' The energy value used in the calculations (in the same units as the Hamiltonian).
988 %> @
param 'EF' The Fermi energy in the same units as the Hamiltonian. (overrides the one comming from the external source)
989 %> @
param 'silent' Set true for suppress the output messages.
990 %> @
param 'transversepotential' A function
handle pot=f(
#coordinates coords) to calculate the transverse potential in the cross section of the ribbon. 991 %> @
param 'leadmodel' A
function handle #
Lead clead=f( idx, E, varargin ) of the alternative lead model with equivalent inputs and
return values as
class #
Lead and with
E standing for the energy.
992 %> @
param 'interfacemodel' A
function handle f( #
InterfaceRegion ) to manually adjus the
interface regions. (Usefull when 'leadmodel' is also given.)
993 %> @
param 'Opt' An instance of the structure #
Opt. (Overrides data in the input file
if given)
994 %> @
param 'param' An instance of the structure #
param. (Overrides data in the input file
if given)
995 %> @
param 'q' The transverse momentum quantum number.
996 %> @
param 'B' The strength of the magnetic field in the units of Tesla.
997 %> @
param 'hole' An instance of structure #
hole.
999 %> @
param 'radius' The radius of the
hole in units of lattice constant.
1000 %> @
param 'rCC' The atomic distance.
1001 function InputParsing(obj, varargin)
1004 p.addParameter(
'width', obj.width);
1005 p.addParameter(
'height', obj.height);
1006 p.addParameter(
'filenameIn', [pwd, filesep,
'Basic_Input_zigzag_leads.xml']);
1007 p.addParameter(
'filenameOut', [pwd, filesep,
'Basic_Input_running_parameters.xml']);
1008 p.addParameter(
'WorkingDir', pwd);
1009 p.addParameter(
'E', obj.E);
1010 p.addParameter(
'EF', obj.EF);
1011 p.addParameter(
'silent', obj.silent);
1012 p.addParameter(
'transversepotential', obj.transversepotential);
1013 p.addParameter(
'leadmodel', obj.leadmodel); %individual physical model
for the contacts
1014 p.addParameter(
'interfacemodel', obj.interfacemodel); %individual physical model
for the
interface regions
1015 p.addParameter('
Opt', []);
1016 p.addParameter(
'param', []);
1017 p.addParameter(
'q', []);
1019 p.addParameter(
'B', obj.B);
1020 p.addParameter(
'hole', obj.hole);
1021 p.addParameter(
'scatterers', obj.scatterers);
1022 p.addParameter(
'radius', 0);
1023 p.addParameter(
'rCC', 1.42);
1025 p.parse(varargin{:});
1027 InputParsing@
Ribbon( obj,
'width', p.Results.width, ...
1028 'height', p.Results.height, ...
1029 'filenameIn', p.Results.filenameIn, ...
1030 'filenameOut', p.Results.filenameOut, ...
1031 'WorkingDir', p.Results.WorkingDir, ...
1032 'E', p.Results.E, ...
1033 'EF', p.Results.EF, ...
1034 'silent', p.Results.silent, ...
1035 'transversepotential', p.Results.transversepotential, ...
1036 'leadmodel', p.Results.leadmodel, ...
1037 'interfacemodel', p.Results.interfacemodel, ...
1038 'Opt', p.Results.Opt, ...
1039 'param', p.Results.param, ...
1043 obj.B = p.Results.B;
1044 obj.hole = p.Results.hole;
1045 obj.scatterers = p.Results.scatterers;
1046 obj.rCC = p.Results.rCC;
1048 if p.Results.radius > 0
1049 obj.hole.radius = p.Results.radius;
1052 obj.waveFncDirnameFull = [pwd, filesep,
'wave_function_full'];
1053 obj.waveFncSubDirname = [
'width',num2str(obj.width),
'_height',num2str(obj.height),
'_radius',num2str(obj.hole.radius),
'_B',num2str(obj.B)];
1060 end % methos
protected Structure hole contains the data about the antidot used in class antidot.
siteindexes
A vector of the site indexes.
Structure Opt contains the basic computational parameters used in EQuUs.
A class for calculations on a ribbon of finite width for equilibrium calculations mostly in the zero ...
function CreateHandlesForMagneticField(B_loc, cRibbon_loc)
Creates and set function handles of the magnetic vector potentials in the Ribbon class.
function Transport(Energy, B)
Calculates the conductance at a given energy value.
function CreateRibbon()
gauge transformation on the inverse Green's function
function Hamiltonians(varargin)
Function to create the custom Hamiltonians for the 1D chain.
A class to perform transport calculations on a graphene antidot (i.e., a hollow in a ribbon)....
function createOutput(filename, Opt, param)
This function creates an output file containing the running parameters.
Structure scatterers contains data on the scattering impurities used in class antidot.
A class describing the interface region between the scattering region and a lead.
A class to calculate the Green functions and self energies of a translational invariant lead The nota...
Property a1
A lattice vector in the hexagonal lattice.
function CreateHandlesForMagneticField()
Creates function handles of the vector potentials and apply the magnetic filed in the ribbon Hamilton...
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.
function gauge_field(x, y, eta_B, Aconst, height, lattice_constant)
Scalar gauge field connecting Landaux and Landauy gauges.
Property E
The energy value for which the TrukkosSajatertekek eigenvalue problem was solved.
function createVectorPotential(B_loc)
Creates the function handle of the magnetic vector potential.
function structures(name)