Mailman 3 August 2020 - Kwant-discuss

optimal way to compute LDOS for large sample
by Sergey Slizovskiy 20 Aug '20

20 Aug '20

Dear Colleagues, (sorry, this is a repost of the same question to a new thread) I have a taks to compute local density of states at a single point for a large homogeneous sample with narrow STS tip potential and in magnetic field. I have tried several methods in Kwant to do this: 1) Naive diagonalization of H for system with no leads and direct calculation of LDOS. It works, but not for large systems, so, I see strong finite-size effects. 2) Sparse Lanczos method with several eigenvalues/functions extracted near each point in the grid of energy values. 3) KPM method. I am not happy with the energy resolution I can get with it. 4) Take a smaller main sample (with the size of the tip potential), but avoid the finite-size effects by attaching a lot of leads pointing in , e.g. 6 directions. Then, use the kwant ldos function to scan LDOS over a grid of energies. Use magnetic_gauge to implement B field in both sample and scattering region. I see the last possibility as the quickest, but it gives me something that seems wrong to me. Does kwant ldos function work correctly if all the leads are insulating for a given energy? Does it correctly catch the localized states in the scattering region? Thank you for any feedback, Sergey Below is a copy of a reply by Pablo Piskunow Dear Sergey, In terms of scaling, the most convenient is the KPM method, since it scales linearly with the inverse energy resolution (that is, the number of moments), and linearly with the system size. Could you clarify what is the issue with the approach 3)? I am not sure what is the energy resolution that you want to achieve. You can conveniently set this value when creating a `kwant.kpm.SpectralDensity` instance via the arguments `energy_resolution` or (exclusive) `num_moments`. The latter is related to the energy resoution by `\delta \pi / N`, where `\delta` is the full width of the spectrum, and `N` the size of your system. In the case of a 2D, that will be similar to `L^2`, where `L` is the linear size. Furthermore, you could progressively increase the energy resolution by adding moments: ``` spectrum = kwant.kpm.SpectralDensity(...) spectrum.add_moments(...) ``` If you want to resolve individual energy levels then the approach 1) and 3) will have the same scaling, however, the KPM method will not give you eigenvalues. Finally, since you want to get the LDOS at a single site or a small set of sites near the STS tip, you should use the `kwant.kpm.LocalVectors` vector factory, setting the `where` argument to the desired spot. ``` def center(site): return site.tag == (0, 0) vector_factory = kwant.kpm.LocalVectors(syst, where=center) spectrum = kwant.kpm.SpectralDensity(syst, num_moments=100, num_vectors=None, vector_factory=vector_factory) ``` Note that `num_vectors=None` will let the `kwant.kpm` module to figure out how many vectors does the vector factory produce. Otherwise, it should be at most the*total number of orbitals* defined by the `where` function, that is, sites times orbitals per site. As you can see, the `kwant.kpm` module is the most suited for this problem, specially when the sample is very large. I hope this helps. Best regards, Pablo n the meantime, you can create a `spectrum` instance, and evaluate the KPM expansion at any energy, for example by ``` spectrum = kwant.kpm.SpectralDensity(...) energy_array = np.linspace(0, 1, 200) density_array = spectrum(energy_array) ``` Another note, if the Hamiltonian is huge, you will save some overhead by passing the bounds of the spectrum explicitly: `bounds=(e_min, e_max)`. These values must be strictly below and above the edges of the spectrum, otherwise the KPM expansion diverges. ah, I forgot. The automatic energy points (`spectrum.energies`) are fixed once you fix the bounds of the spectrum. So if for a range of magnetic fields you pass the same bounds, then the energies will be at the same points and so the densities will be easier to plot or compare. The width of the spectrum might change for different values of the magnetic field; and even if not, there is small randomness associated with finding the bounds, that you can get rid of by passing a seed to the random number generator like `rng=0`. If you choose to provide the bounds explicitly, you can disregard this random number generator stuff.

2 1

Andreev reflection in Graphene
by Martin Fonnum Jakobsen 20 Aug '20

20 Aug '20

Hello, I have recently started using KWANT. I hope this is the correct way to ask questions. I'm interested in reproducing Fig. 4 in PRL 97, 067007 (2006) - Specular Andreev Reflection in Graphene. This has been done in KWANT in the PhD thesis of Tibor Sekera - "Quantum transport of fermions in honeycomb lattices and cold atomic systems", chapter 4. The problem is similar to tutorial 2.6. In my problem, I have a tight-binding Hamiltonian living on a hexagonal lattice, from which I can construct a discrete BdG equation (See attachment). In addition, to electron-hole degrees of freedom I'm also interested in including spin because later I will add magnetic fields and such. I consider a normal metal lead attached to a superconducing scattering region with zero barrier at the interface. Since I have both electron-hole and spin degrees of freedom my BdG equation is 4x4. To account for electron-hole degrees of freedom I have specified that at each site there should be two orbitals. In the tutorial 2.6 it is emphasized that it is important to specify the correct conservation law in the N lead when we want to calculate the conductance. My question is how should the conservation law (and particle hole symmetry) in the normal-metal lead look like for my case? I have tried the standard expression L_0 = kwant.Builder(sym,conservation_law=s_z) where s_z is the Pauli matrix [[1,0],[0,-1]]. In this case I obtain a dimension mismatch error. I have also tried a generalized 4x4 matrix version of this L_0 = kwant.Builder(sym,conservation_law=tinyarray.array([[1,0,0,0],[0,-1,0,0],[0,0,1,0],[0,0,0,-1]])) Or L_0 = kwant.Builder(sym,conservation_law=tinyarray.array([[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]])). In these cases I obtain the error "Single `onsite` matrix of shape (4, 4) provided but there are 2 orbitals per site in the system". I think my issue is that I do not fully understand what the Conservation law actually does. Initially I thought the conservation law was a matrix C such that C.H.C^-1 should return a block diagonal matrix, where the eigenvalues of C refers to each block. Is this not an accurate interpretation? I have attached a PDF with the BdG equation, and my python code as a .py and .txt file. The content of the .py and .txt file are identical, which format do you prefer in the future? I hope this is OK, and not too overwhelming. I also tried to subscribe (with this e-mail) to the mailing list, but I'm not sure that it worked. Yours sincerely, Martin F Jakobsen

5 6

Is Kwant capable of systems that have different site with different orbitals per site.
by Xiangru Kong 20 Aug '20

20 Aug '20

I know there is Tbmodels package can deal with wannier90 output tight binding files. But I think TBmodels group all the orbitals into one site. I was trying to do things that can assign hoppings and on-site energies in "Real" systems that have differnet sites but with different orbitals per site. The errors are the following when I try to use kwant_model.hamiltonian_submatrix module --------------------------------------------------------------------------- ValueError Traceback (most recent call last) /Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/builder.py in hamiltonian(self, i, j, params, *args) 1885 try: -> 1886 value = value(site, *args) 1887 except Exception as exc: /Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/wraparound.py in f(*in_args) 154 if callable(val): --> 155 acc = acc + val(*out_args) 156 else: ValueError: operands could not be broadcast together with shapes (6,6) (10,10) The above exception was the direct cause of the following exception: UserCodeError Traceback (most recent call last) <ipython-input-34-1b24fb1489bf> in <module>() 3 params={key: val for key, val in zip(['k_x', 'k_y', 'k_z'], 2 * np.pi * np.array(k_list[0]))} 4 ) ----> 5 ) for k in k_list] <ipython-input-34-1b24fb1489bf> in <listcomp>(.0) 3 params={key: val for key, val in zip(['k_x', 'k_y', 'k_z'], 2 * np.pi * np.array(k_list[0]))} 4 ) ----> 5 ) for k in k_list] /Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/_common.py in inner(*args, **kwargs) 70 if sig.bind(*args, **kwargs).arguments.get(parameter_name): 71 warn() ---> 72 return f(*args, **kwargs) 73 74 return inner kwant/_system.pyx in kwant._system.hamiltonian_submatrix() /Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/builder.py in hamiltonian(self, i, j, params, *args) 1892 ', '.join(map('"{}"'.format, missing))) 1893 raise TypeError(''.join(msg)) -> 1894 _raise_user_error(exc, value) 1895 else: 1896 edge_id = self.graph.first_edge_id(i, j) /Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/builder.py in _raise_user_error(exc, func) 1814 msg = ('Error occurred in user-supplied value function "{0}".\n' 1815 'See the upper part of the above backtrace for more information.') -> 1816 raise UserCodeError(msg.format(func.__name__)) from exc 1817 1818 UserCodeError: Error occurred in user-supplied value function "f". See the upper part of the above backtrace for more information.

2 3

Unit of current density in current plot
by Prof. CHAN Kwok Sum 19 Aug '20

19 Aug '20

Dear Kwant Developers, I would like to know the unit of the current density shown in the current plot. Does it always equal "unit of current/relwidth" for a 2D system? I understand that abswidth can replace relwidth. If I do not specify any abswidth, will the relwidth be replaced? Finally, the relwidth is defined as , as a fraction of the length of the longest side of the bounding box. What is the definition of the bounding box? Is it the area in which the current is plotted? Thank you in advance for your help. Regards, KS Chan Disclaimer: This email (including any attachments) is for the use of the intended recipient only and may contain confidential information and/or copyright material. If you are not the intended recipient, please notify the sender immediately and delete this email and all copies from your system. Any unauthorized use, disclosure, reproduction, copying, distribution, or other form of unauthorized dissemination of the contents is expressly prohibited.

2 1

Hall effect with real size graphene
by Khani Hosein 14 Aug '20

14 Aug '20

Dear all, I am trying a graphene bulk with magnetic field, and I use spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None, vector_factory=s_factory, rng=0) energies, densities = spectrum() to calculate the dos and compare with some previous studies. My definition of the hopping is: def Hop_magnetic(site1,site2): B_magneic=200. # units Tesla a0=0.142 # nm x1,y1=site1.pos x2,y2=site2.pos xy=0.5*(x1+x2)*sqrt(3)*a0*(y1-y2)*sqrt(3)*a0 # the units is nm**2 phb=1j*2.*pi*B_magneic/4.135667 # phi0=h/e=4.135667*e-15 V*s return t*exp(xy*phb*0.001) I found that the effect of the magnetic field is too small and I can not obtain the correct results. Could you give me some suggestions? Best regards Khani My code is pasted: import kwant from matplotlib import pyplot from numpy import sqrt,pi,exp import numpy as np def make_system(r=8, t=1, tp=-0.1): lat = kwant.lattice.honeycomb(norbs=1) a, b = lat.sublattices def circle(pos): x, y = pos return x**2 + y**2 < 100**2 #-100<x<100 and -100<y<100 def Hop_magnetic(site1,site2): B_magneic=200. # the unit of B_magneic is Tesla=V*s/m**2 a0=0.142 # nm x1,y1=site1.pos x2,y2=site2.pos xy=0.5*(x1+x2)*sqrt(3)*a0*(y1-y2)*sqrt(3)*a0 # the units change it to nm**2 phb=1j*2.*pi*B_magneic/4.135667 # phi0=h/e=4.135667*e-15 V*s, the unit of B_magneic is Tesla=V*s/m**2 return t*exp(xy*phb*0.001) def Nearest_Hop(site1,site2): return Hop_magnetic(site1,site2) syst = kwant.Builder() syst[lat.shape(circle, (0, 0))] = 0 syst[lat.neighbors()] = Nearest_Hop syst.eradicate_dangling() return lat, syst.finalized() lat, fsyst = make_system() #kwant.plot(fsyst) where = lambda s: np.linalg.norm(s.pos) < 1 s_factory = kwant.kpm.LocalVectors(fsyst, where) spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None, vector_factory=s_factory, rng=0) energies, densities = spectrum() pyplot.figure() pyplot.plot(energies, densities) pyplot.xlabel("energy [t]") pyplot.ylabel("DOS") pyplot.show()

2 1

Width of leads larger than the width of the scattering region
by octaviors＠gmail.com 13 Aug '20

13 Aug '20

I am using two leads (left and right) where their widths are larger the scattering region. For instance, the sites of the scattering region are between the vertical interval y1 <= y <= y2, so its width is W = y2-y1 +1. The sites of the leads are of course between 0 <= y <= W_lead. So, W < W_lead. The problem is that attach_lead doesn't work, unless one includes sites in the scattering region for 0<= y < y1, y2<y<=W_lead. Thus I did not define hoppings for those zones, however attach_lead() imposes spurious horizontal hoppings at the frontiers, otherwise it gives an error message. it seems that attach_lead needs to connect the whole frontier of any lead. Is there any way for avoiding these spurious horizontal hoppings ?.

2 1

discretize continuous Hamiltonians
by Khani Hosein 08 Aug '20

08 Aug '20

Dear authors, I want to discretize continuous Hamiltonians of Weyl semimetal (PRL 115,246603 (2015)) using Kwant, but we have sin(kz),sin(ky)...in the Hamiltonian, so it always shows the error message. I have already imported sin and cos from math. Could you give me any suggestion? Hosein Khani

4 8

Migrating spin_orbit.py to real unit system
by Nathaniel Gregory Seil 06 Aug '20

06 Aug '20

Hello, I would like to model a 2DEG under a magnetic field and have been investigating the spin_orbit.py code. I would like to modify this code so that it utilises real units and models upon a to-scale lattice. However, I am a little bit unclear as to how to progress. I noticed in the base code we have a = 1, t = 1 and hbar = 2 (seen from the Pauli spin matrices which normally have hbar/2 as a coefficient in front). To represent the real life scenario of a lattice with constant a = 5E-10 m, would we do the following: * a = 5E-10 m * hbar = 1.05457E-34 J s (standard SI units) * m_electron = 9.109E-31 kg (standard SI units) * t = hbar ^ 2 / (2 * m_electron * a^2) and then just substitute these values into the rest of the code? I get the result of 0 conductivity for the lead voltages -0.046 to 0.105V (code shows "energies=[0.01 * i*t - 0.3*t for i in range(100)]"). I also get this for -0.0046 to 0.0105V and for -0.46 - 1.05V. Do these results seem right? I have attached my code with these substitutions and the results. Many thanks, Nathaniel Seil

2 5

Re: Peierls phase from kwant.physics.magnetic_gauge
by Pawel Wojcik 02 Aug '20

02 Aug '20

Dear Pablo, Yes, now, everything works. Thanks a lot for Your help. Best Regards, Pawel >>> <pablo.perez.piskunow(a)gmail.com> 07/31/20 10:13 PM >>> Hi again Paweł, I realized that `substutited` *returns* a new builder, and does not change it in-place. The docs of `magnetic_gauge` were misleading there. You can simply replace ``` l_lead = l_lead.substituted(peierls='peierls_llead') ... r_lead = r_lead.substituted(peierls='peierls_llead') ``` in your code, and should work.

1 0