1 %%
Transport calculation through a graphene QD
using self-consistent potential - based on EQuUs v4.8
2 % Copyright (C) 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 Examples
20 %> @brief Example to calculate transport through a graphene quantum dot
using self-consistent potential. Compare results to Fig.3 in 2014 Nanotechnology 25 465201. In these calculations we use a leads defined on a square lattice and
do the
self-consistent potential iterations only
for \f$E_F =0\f$.
21 %> @
param filenum The identification number of the filenema
for the exported data (
default is 1).
23 %> S. Krompiewski, Nanotechnology 25 465201 (2014).
25 %> EQuUs v4.9 or later
31 %> @brief Example to calculate transport through a graphene quantum dot
using self-consistent potential. Compare results to Fig.3 in 2014 Nanotechnology 25 465201. In these calculations we use a leads defined on a square lattice and
do the
self-consistent potential iterations only for \f$E_F =0\f$.
32 %> @
param filenum The identification number of the filenema
for the exported data (
default is 1).
35 if ~exist(
'filenum',
'var')
39 filename = mfilename('fullpath');
40 [directory, fncname] = fileparts( filename );
42 % filename containing the input XML
43 inputXML = 'Input_Transport.xml';
44 % Parsing the input file and creating data
structures 46 vargamma =
param.scatter.vargamma;
48 % Strength of the self-consitent intercation in the Stoner model
50 % Zeeman energy for the initial iteration to bring the system out of the trivial nonmagnetic solution
51 Zeeman_energy = 0.1*vargamma;
54 % number of energy point on the contour
58 % width of the junction
60 % length of the junction
62 % filename containing the output XML
63 outfilename = [fncname,'_',num2str( filenum )];
66 % false to supress output messages
68 % Using self energy in the transport calculations
72 % noninteracting transport results
73 Evec_non_interacting = [];
74 Conductance_non_interacting = [];
76 % self_consistent transport results
77 Conductance_self_consistent = [];
78 Evec_self_consistent = [];
80 % Class for self_consistent iterations
83 % creating output directories
86 % Calculates the non-interacting transport through a graphene island
89 % Calculates the self-consistent transport through a graphene island
95 % export the calculetd data
96 save( [outputdir,'/',outfilename,'.mat']);
99 %> @brief calculates the transport with self-consistent potential
102 display('****** Calculating self-consistent transport ******')
105 % Creates the class for self-consistent iterations for magnetization up
106 cSelfConsistent =
SelfConsistent(
Opt, 'magnetizaion_direction', [], 'tolerance', 2e-3, 'newIterationWeight', 0.3, 'outFileName', '
test.mat', 'MaxIter', 500, 'plotIterations', false, 'resume', false ) ;
111 Evec_self_consistent = Emin:(Emax-Emin)/Enum:Emax;
112 Conductance_self_consistent = zeros(size(Evec_self_consistent));
114 % calcualting the self-consistent potential
116 % opening the parallel pool
119 poolmanager.openPool();
121 % running the self-consistent iterations
124 % closing the parallel pool
125 poolmanager.closePool();
127 for idx = 1:length(Evec_self_consistent)
129 Conductance_self_consistent(idx) =
Transport(Evec_self_consistent(idx), true);
131 display(['Energy = ', num2str(Evec_self_consistent(idx)), ', conductance = ', num2str(Conductance_self_consistent(idx))]);
141 %> @brief calculates the non-interacting transport
144 display('****** Calculating transport in the non-interacting case ******')
147 % oreallocating arrays for the results
151 Evec_non_interacting = Emin:(Emax-Emin)/Enum:Emax;
152 Conductance_non_interacting = zeros(size(Evec_self_consistent));
154 for idx = 1:length(Evec_non_interacting)
155 Conductance_non_interacting(idx) =
Transport(Evec_non_interacting(idx), false);
156 display(['Energy = ', num2str(Evec_non_interacting(idx)), ', conductance = ', num2str(Conductance_non_interacting(idx))])
167 %> @brief Calculates the conductance for a given energy for both self-consistent and non-interacting cases.
168 %> @
param Energy The enrgy value
169 %> @
param self_consistent Logical value. Set true to use the self-consistent potential, or false to do the non-interacting calculations.
170 %> @return Returns with the calculated conductance in units of \f$e^2/\hbar \f$.
175 % creating the
Ribbon interface for the SNS structure
176 cRibbon =
Ribbon('width', width, 'height', height, 'EF', EF, 'silent', silent, '
Opt',
Opt, '
param',
param, 'filenameOut', fullfile(outputdir, [outfilename, '.xml']), ...
179 % setting the energy value
180 cRibbon.setEnergy( Energy )
182 % calculating the Green
function of the scattering reigon
183 cRibbon.CalcFiniteGreensFunctionFromHamiltonian('PotInScatter', @(CreateH,E)
Pot_Self_Consistent(CreateH, E, self_consistent), 'onlyGinv', true);
185 % calculating the leads
186 coordinates_shift = [-2, height+1];
187 recalculateLeads = [1 2];
188 cRibbon.FL_handles.LeadCalc('coordinates_shift', coordinates_shift, 'SelfEnergy', useselfEnergy, 'SurfaceGreensFunction', ~useselfEnergy, 'leads', recalculateLeads, 'leadmodel', @
Square_leads);
190 % creating interface regions
191 for idx = 1:length(cRibbon.Interface_Regions)
192 cRibbon.CreateInterface( idx );
197 Dysonfunc = @()(cRibbon.CustomDysonFunc( 'recalculateSurface', [], 'constant_channels', true, 'decimate', true, 'SelfEnergy', useselfEnergy ));
198 cRibbon.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
200 % calculating the conductance
201 ret = cRibbon.FL_handles.Conductance2();
208 %> @brief Calculates the self-consistent potential
210 % reset the self-consistent class to forget previous iterations, but keep the last calculated density
211 cSelfConsistent.Reset();
213 % Do the self-consistent iterations
214 cSelfConsistent.runSelfConsistentIterations( @
DensityCalc );
220 %> @brief Calculates the density using the predefined graphene lattice framework (
Ribbon class)
221 %> @return [1] The calculated density data
222 %> @return [2] An instance of structure
#junction_sites describing the sites in the calculated density 225 % creating the
Ribbon interface for the SNS structure
226 cRibbon =
Ribbon(
'width', width,
'height', height,
'EF', EF,
'silent', silent,
'Opt',
Opt,
'param',
param,
'filenameOut', fullfile(outputdir, [outfilename,
'.xml']), ...
230 cDensity =
Density(
Opt,
'T', T,
'scatterPotential', @(CreateH, E)(
Pot_Self_Consistent(CreateH, E,
true)),
'junction', cRibbon,
'useSelfEnergy', useselfEnergy);
232 % Calculate the density
233 [density_data,
junction_sites] = cDensity.DensityCalc(
'Edb', Edb,
'Emin', -11.0286 );
239 %> @brief Create a potential from a
self-consistent data, and removing unecessary
sites from the scattering region
241 %> @
param E The energy value
242 %> @
param SelfConsistent Logical value. Set
true to calculate the
self-consistent potential, or
false otherwise
249 % determine
sites to be removed to reproduce the scattering region in 2014 Nanotechnology 25 465201
250 indexes2remove = ~(y<max_y);
252 CreateH.Write(
'kulso_szabfokok', []);
254 y( indexes2remove ) = max_y/2;
259 kulso_szabfokok_logical = ~(y<max_y & y>min_y);
265 % removing
sites from the scattering region
266 CreateH.RemoveSites( indexes2remove )
277 % adding the
self-consistent interaction
281 if isempty( cSelfConsistent.density)
284 % adding Zeeman term at the edges the first iteration to bring the
285 %system put of the trivial solution
290 %ret = (exp(-(x-min_x).^2/(2*lattice_constant^2)) - exp(-(x-max_x).^2/(2*lattice_constant^2)) )*Zeeman_energy;
291 ret = -ones(size(x))*Zeeman_energy;
296 % Stoner interaction model of Nanotechnology 25 (2014) 465201
297 scatter_density = cSelfConsistent.density(cSelfConsistent.
junction_sites.Scatter.site_indexes);
302 ret(
coordinates.spinup) = 0.5*U*(density_scatter_spinup-density_scatter_spindown);
303 ret(~
coordinates.spinup) = -0.5*U*(density_scatter_spinup-density_scatter_spindown);
310 %> @brief Create square leads
Hamiltonians with the same inputs and output as
#Transport_Interface.SurfaceGreenFunctionCalculator. 312 p_inner = inputParser;
314 p_inner.addParameter(
'createCore', 0);
315 p_inner.addParameter(
'Just_Create_Hamiltonians', 0);
316 p_inner.addParameter(
'shiftLead', 0);
317 p_inner.addParameter(
'coordinates_shift', 0);
318 p_inner.addParameter(
'transversepotential', []);
319 p_inner.addParameter(
'Lead',[]);
320 p_inner.addParameter(
'gauge_field', [] );% gauge field
for performing gauge transformation
321 p_inner.addParameter(
'SelfEnergy',
false); % set
true to calculate the
self energy of the semi-infinite lead
322 p_inner.addParameter(
'SurfaceGreensFunction',
true );% set
true to calculate the surface Greens
function of the semi-infinite lead
323 p_inner.addParameter(
'q', [] ) %transverse momentum
327 createCore = p_inner.Results.createCore;
328 Just_Create_Hamiltonians = p_inner.Results.Just_Create_Hamiltonians;
329 shiftLead = p_inner.Results.shiftLead;
330 coordinates_shift = p_inner.Results.coordinates_shift;
331 transversepotential = p_inner.Results.transversepotential;
332 Lead = [];%p_inner.Results.Surface_tmp;
333 gauge_field = p_inner.Results.gauge_field; % gauge field
for performing gauge transformation
334 SelfEnergy = p_inner.Results.SelfEnergy;
335 SurfaceGreensFunction = p_inner.Results.SurfaceGreensFunction;
336 q = p_inner.Results.q;
339 Opt2.Lattice_Type =
'Square';
344 param2.Leads{idx}.epsilon =
param.
Leads{idx}.epsilon;
345 param2.Leads{idx}.vargamma =
param.
Leads{idx}.vargamma;
346 param2.Leads{idx}.M =
width;
351 Lead_ret = FL_handles.SurfaceGreenFunctionCalculator(idx,
'createCore', createCore, ...
352 'Just_Create_Hamiltonians', Just_Create_Hamiltonians, ...
353 'shiftLead', shiftLead, ...
354 'coordinates_shift', coordinates_shift, ...
355 'transversepotential', transversepotential, ...
357 'SelfEnergy', SelfEnergy, ...
358 'SurfaceGreensFunction', SurfaceGreensFunction, ...
365 %> @brief Method to adjust the
Hamiltonians of the
interface regions
366 %> @
param Interface_Region An instance of
class #Interface_Region.
369 Kcoupling = Interface_Region.Read(
'Kcoupling');
370 Kcouplingadj = Interface_Region.Read(
'Kcouplingadj');
373 [rows, cols] = find( Kcoupling );
376 Kcoupling = Kcoupling(rows, :);
377 Kcouplingadj = Kcouplingadj(:, rows);
379 Interface_Region.Write(
'Kcoupling', Kcoupling);
380 Interface_Region.Write(
'Kcouplingadj', Kcouplingadj);
389 % creating figure in units of pixels
390 figure1 = figure(
'Units',
'Pixels',
'Visible',
'off');
392 % font size on the figure will be 16 points
395 % setting the limits of the axis
398 % creating axes of the plot
399 axes1 = axes(
'Parent',figure1, ...
401 'FontSize', fontsize,...
404 'Units',
'Pixels', ...
405 'FontName',
'Times New Roman');
409 xlabel(
'E [eV]',
'FontSize', fontsize,
'FontName',
'Times New Roman',
'Parent', axes1);
412 ylabel(
'G [e^2/h]',
'FontSize', fontsize,
'FontName',
'Times New Roman',
'Parent', axes1);
414 % reserve a variable
for legend titles
417 % plot the non-interacting and
self-consistent data
418 if ~isempty( Conductance_non_interacting )
419 plot(Evec_non_interacting, real(Conductance_non_interacting),
'b',
'LineWidth', 2,
'Parent', axes1)
420 legend_labels{end+1} =
'U = 0';
423 if ~isempty( Conductance_self_consistent )
424 plot(Evec_self_consistent, real(Conductance_self_consistent),
'rx',
'Parent', axes1)
425 plot(Evec_self_consistent, real(Conductance_self_consistent),
'r',
'LineWidth', 1.5,
'Parent', axes1)
426 legend_labels{end+1} = [
'U = ', num2str(U/vargamma),
'\gamma'];
431 % plot the non-interacting and
self-consistent dat
432 fig_legend = legend(axes1, legend_labels);
433 set(fig_legend,
'FontSize', fontsize*0.8,
'FontName',
'Times New Roman',
'Box',
'off',
'Location',
'NorthEast')
435 % set the paper position in releases greater than 2016
436 ver = version(
'-release');
437 if str2num(ver(1:4)) >= 2016
438 % setting the position and margins of the plot, removing white spaces
439 figure1.PaperPositionMode =
'auto';
440 fig_pos = figure1.PaperPosition;
441 figure1.PaperSize = [fig_pos(3) fig_pos(4)];
444 set(axes1, 'Position', get(axes1, 'OuterPosition') - get(axes1, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1]);
445 Position_figure = get(axes1, 'OuterPosition');
446 set(figure1, 'Position', Position_figure);
450 print('-depsc2', fullfile(outputdir,[outfilename, '.eps']))
451 print('-dpdf', fullfile(outputdir,[outfilename, '.pdf']))
462 % sets the output directory
464 resultsdir = 'results';
465 mkdir(fullfile(directory, resultsdir) );
466 outputdir = resultsdir;
function test(arg1, arg2)
Brief description of the function.
lead_param Leads
A list of structures lead_param containing the physical parameters for the scattering region.
function Transport(Energy, self_consistent)
Calculates the conductance for a given energy for both self-consistent and non-interacting cases.
Property width
The number of sites in the cross section.
A class for controlling the parallel pool for paralell computations.
Structure Opt contains the basic computational parameters used in EQuUs.
function Square_leads(idx, E, varargin)
Create square leads Hamiltonians with the same inputs and output as Transport_Interface....
A class for calculations on a ribbon of finite width for equilibrium calculations mostly in the zero ...
Property coordinates
An instance of the structure coordinates.
function CalcSelfConsistentTransport()
calculates the transport with self-consistent potential
A class to evaluate self-consistent potentials depending on electron (spin resolved) densities.
A class to evaluate the Dyson equation and to calculate the scattering matrix for equilibrium process...
function Hamiltonians(varargin)
Function to create the custom Hamiltonians for the 1D chain.
Property q
A transverse momentum.
function DensityCalc()
Calculates the density using the predefined graphene lattice framework (Ribbon class)
Property kulso_szabfokok
A list of the sites to be kept after decimation.
function CalcSelfconsistentPotential()
Calculates the self-consistent potential.
A class to calculate the Green functions and self energies of a translational invariant lead The nota...
Property varargin
list of optional parameters (see http://www.mathworks.com/help/matlab/ref/varargin....
A class to calculate the onsite desnity useful for self-consistent calculations.
Class containing physical parameters of a particular lead defined on a square lattice.
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 InterfaceModel(Interface_Region)
Method to adjust the Hamiltonians of the interface regions.
Property Opt
An instance of structure Opt.
function CalcTransport()
calculates the non-interacting transport
Property param
An instance of structure param.
function gauge_field(x, y, eta_B, Aconst, height, lattice_constant)
Scalar gauge field connecting Landaux and Landauy gauges.
function Read(MemberName)
Query for the value of an attribute in the class.
function Self_Consistent_Transport(filenum)
Example to calculate transport through a graphene quantum dot using self-consistent potential.
function structures(name)
function Pot_Self_Consistent(CreateH, E, SelfConsistent)
Create a potential from a self-consistent data, and removing unecessary sites from the scattering reg...
Structure junction_sites contains data to identify the individual sites of the leads,...
A class to create and store Hamiltonian of the scattering region.