Re: [Kwant] gas adsorption by a crystal

[Anton Akhmerov]

Hi Marco,

I'm not familiar with Kubas bonding model, can you please provide more
details and background about what you have in mind?

Best,
Anton

On Sat, 18 Apr 2020 at 11:51, Marco Di Gennaro (TME)
 wrote:
>
> Dear everyone,
>
>
>
> I would like to study gas dynamics and absorption within a transition-metal 
> based crystal.
>
> Is possible to realize a simplified Kubas bonding model in tight binding 
> through a simple s-d interaction?
>
> Could I distinguish physisorption by chemisorption by Kubas adsorption at 
> this level of theory?
>
>
>
> Any comment or suggestion would be highly appreciated.
>
> Thank you
>
>
>
> -
>
> Marco Di Gennaro
>
> Toyota Motor Europe
>
> -
>
>
>
> This e-mail may contain confidential information. If you are not an addressee 
> or otherwise authorised to receive this message, you should not use, copy, 
> disclose or take any action based on this e-mail. If you have received this 
> e-mail in error, please inform the sender promptly and delete this message 
> and any attachments immediately.

Re: [Kwant] Density of states for a lead

[Anton Akhmerov]

Dear Jacopo,

There's nothing out of the box, however the calculation is rather
straightforward. If you want to compute the LDOS of a lead at a given
energy, compute its modes (via lead.modes), then take the propagating
ones and square their sum.
Since the lead modes are normalized to carry unit current, this
procedure is sufficient. If you only wan the overall DOS, sum the
inverse velocities, as returned by the modes. Come to think of it, it
would make a lot of sense to implement density of states as a
convenience method of propagating modes. I've opened
https://gitlab.kwant-project.org/kwant/kwant/-/issues/376 ; anyone
interested?

Best,
Anton

On Thu, 16 Apr 2020 at 22:07, Jacopo Settino  wrote:
>
> Dear all,
> is it implemented a way to calculate the density or the density of states of 
> a lead (defined as kwant.TranslationalSymmetry ...)?
> Thank you.
> Regards,
> Jacopo
>
>

Re: [Kwant] [EXTERNAL] Re: Stability of the retarded Green's function calculation

[Anton Akhmerov]

Hey Joe,

> When you say "computing a self-energy corresponds to partially fixing the 
> order of Gaussian elimination" do you mean that computing the self-energy 
> forces us to invert Φ, and *this* is what leads to an instability?
> As I write this I realise that I am just restating what you said, but I would 
> like a bit of clarity on this point.

Indeed: when there's a lead bound state at the energy of interest, Φ
becomes not invertible. More generally, there are no guarantees on the
condition number of Φ.

Cheers,
Anton

> Thanks,
>
> Joe
>
>
>
> -----Original Message-
> From: Anton Akhmerov 
> Sent: Sunday, March 22, 2020 8:41 AM
> To: Joseph Weston (Aquent LLC - Canada) 
> Cc: kwant-discuss@kwant-project.org
> Subject: [EXTERNAL] Re: [Kwant] Stability of the retarded Green's function 
> calculation
>
> Hi Joe,
>
> Computing the Green's function or the self-energy of the lead fails when a 
> terminated lead has a bound state at the energy we're looking at (so there's 
> an exact pole). For example a lead made out of a topological superconductor 
> always has this problem at zero energy. At the same time the mode solvers 
> face no singularities in this case. In linear algebra terms computing a 
> self-energy corresponds to partially fixing the order of Gaussian 
> elimination, while the most general case may require pivoting incompatible 
> with this order.
>
> That's what concerns stability. The speed difference is, I think,
> insignificant: the number of right hand sides is smaller with the modes 
> solver since we use one per incoming mode, not one per degree of freedom. 
> However in practice obtaining the LU decomposition causes the largest 
> slowdown.
>
> Cheers,
> Anton
>
>
> On Wed, 18 Mar 2020 at 22:12, Joseph Weston (Aquent LLC - Canada) 
>  wrote:
> >
> > Hello again,
> >
> >
> >
> > I Noticed a couple of minor mistakes in my previous email:
> >
> >
> >
> > I used the term “out of the box” linear solvers, when I meant “black box” 
> > linear solvers. By this I meant that we are not tailoring the algorithm 
> > used for solving the linear system based on the properties of the LHS (e.g. 
> > as we do in RGF by the choice of scattering region slices), but rather 
> > trusting that the applied mathematicians have done their job properly, and 
> > any incidental structure in the LHS is found “automatically” by the solver.
> > I claimed that the RHSs for the scattering problem were “indicator vectors 
> > with 1s in [the] extended part”. This is not true, however the 
> > interpretation of the RHS as a “single incoming mode” is correct AFAIK.
> >
> >
> >
> > Happy Kwanting,
> >
> >
> >
> > Joe
> >
> >
> >
> > From: Kwant-discuss  On
> > Behalf Of Joseph Weston (Aquent LLC - Canada)
> > Sent: Wednesday, March 18, 2020 12:08 PM
> > To: kwant-discuss@kwant-project.org
> > Subject: [EXTERNAL] [Kwant] Stability of the retarded Green's function
> > calculation
> >
> >
> >
> > Hello Kwantoptians,
> >
> >
> >
> > I have a question regarding the speed/stability of computing the retarded 
> > Green’s function of a transport setup using Kwant.
> >
> >
> >
> > In the Kwant source-code [1] it is noted that using `kwant.greens_function` 
> > is “often slower and less stable than the scattering matrix calculation”. I 
> > was wondering if you could provide me with some references for this 
> > assertion.
> >
> >
> >
> > I have perused the Kwant paper [2] and (part I of) the thesis of Michael 
> > Wimmer [3] and cannot find any mention of a speed/stability comparison of 
> > the algorithm implemented by `kwant.greens_function` vs `kwant.smatrix`.
> >
> >
> >
> > My understanding is that in both cases a linear system (LHS) is constructed 
> > and solved for different right hand-sides (RHS) using out of the box linear 
> > solvers. The solution for each RHS corresponds to one column of the 
> > retarded Green’s function / extended scattering matrix respectively. The 
> > difference between `kwant.greens_function` and `kwant.smatrix` is then the 
> > following. In the former case the leads are taken into account by added the 
> > retarded self-energy to the LHS and the RHSs are indicator vectors for the 
> > sites of the lead/scattering region interface,  which corresponds to a 
> > “unit impulse” boundary condition. In the latter case the linear system is 
> > “extended” so as to include extra unknowns that correspond to the 
> > scattering amplitud

Re: [Kwant] How to obtain the retarded Green's function for a lead

[Anton Akhmerov]

Hey Joe,

You're absolutely right. One additional reason why a routine for the
lead's Green's function doesn't exist is also that it was not a
priority so far.

It should be possible to obtain the lead's Green's function though. We
would need to take the appropriate linear superpositions of the modes,
which we would then expand (using the identically called function in
modes.py) to the original tight-binding basis.

Cheers,
Anton

On Fri, 13 Mar 2020 at 08:43, Joseph Weston (Aquent LLC - Canada)
 wrote:
>
> Good day,
>
>
>
> I have a question about how best to obtain the retarded Green’s function for 
> a lead in Kwant.
>
>
>
> Kwant already has a way of calculating the retarded self-energy of a lead 
> (the ‘selfenergy’ method of ‘StabilizedModes’ [1]), however there is no 
> equivalent routine for the retarded Green’s function. As far as I am aware 
> one may obtain the self-energy from the Green’s function by sandwiching it 
> between a V (inter-cell hopping) and a V^\dagger, and I was wondering if 
> there is a reason why there is not a routine that returns the Green’s 
> function directly.
>
>
>
> One possibility is because the modes, as calculated by Kwant, are (in 
> general) in a weird basis related to the singular value decomposition of V, 
> and only a part of the SVD is stored with the stabilized modes, which makes 
> it impossible to “get out” the modes in the tight-binding basis as required 
> for the Green’s function.
>
>
>
> Is my understanding correct, or am I missing something?
>
>
>
> Thanks,
>
>
>
> Joe
>
>
>
> [1]: 
> https://gitlab.kwant-project.org/kwant/kwant/-/blob/master/kwant/physics/leads.py#L185

Re: [Kwant] question for energy range in kwant

[Anton Akhmerov]

Hi Daehan,

This happens because the tight binding models have a finite band width.
Compare your results with a band structure of the model you are considering.

Best,
Anton

On Thu, 12 Mar 2020 at 09:00, dae han Park  wrote:

> Hello.
>
> My name is Daehan Park. I am a new user for Kwant.
>
> To accustom me to using Kwant, I start to attack simple quantum mechanics
> problems in kwant program.
>
> In below code, I try to solve a simple square barrier problem.
> ==
> import kwant
> import math
> from matplotlib import pyplot
> syst = kwant.Builder()
>
> wire = kwant.lattice.chain(1)
> t = 1.0   # hopping parameter
> W = 100# the number of atoms in the direction of width
>
> for i in range(W):
> # on-site Hamiltonian
> syst[(wire(i))]=2*t
>
> for i in range(10):
> # on-site Hamiltonian
> syst[(wire(i+10))]=2*t + 0.5
>
> syst[wire.neighbors()] = -t # hopping
>
> # attaching leads
> sys_left_lead = kwant.TranslationalSymmetry([-1])
> left_lead = kwant.Builder(sys_left_lead)
> left_lead[wire(0)] =  2*t
> left_lead[wire.neighbors()] = -t
> syst.attach_lead(left_lead)
> syst.attach_lead(left_lead.reversed())
>
> kwant.plot(syst)
> sys = syst.finalized()
>
> def plot_conductance(sys, energies):
> # Compute conductance
> data = []
> for energy in energies:
> smatrix = kwant.smatrix(sys, energy)
> data.append(smatrix.transmission(1, 0))
>
> pyplot.figure()
> pyplot.plot(energies, data)
> pyplot.xlabel("energy [t]")
> pyplot.ylabel("conductance [e^2/h]")
> pyplot.show()
>
> plot_conductance(sys, energies=[0.01* i + 1e-6 for i in range(500)])
> 
>
> The outcome above the code is as following;
>
> [image: image.png]
>
> My question is why conductance becomes zero when the incident energy is
> over 4.
>
> Thank you for reading my e-mail
>
> Sincerely,
> Daehan Park.
>
> =
> Department of Physics,
> Soongsil University,
> 369, Sangdo-ro, Dongjak-gu,
> Seoul 156-743,
> Republic of Korea (South Korea)
> mail:  m31...@gmail.com
> =
> ᐧ
>

[Kwant] Idea discussion: development internships

[Anton Akhmerov]

Hi everyone!

Together with the other Kwant authors I am planning a grant
application. One activity for which we consider requesting support is
a series of internships—visits to one of our groups with the goal of
developing either a Kwant feature or a related project. These
internships would last for perhaps a couple of months, with the idea
that by the end of the visit, they result in something either
implemented or being close to completion.

We think that this work organization would benefit everyone involved,
and we'm very much interested to explore this idea.

Still, before we commit to that, we would be very much interested to
hear everyone's opinion: Would you or someone you know want to go for
such an internship? What do you think of this idea overall? Any idea
of what project could be useful?

Cheers,
Anton

Re: [Kwant] attach lead

[Anton Akhmerov]

Hi!

Please check the documentation for lattice.shape over here:
https://kwant-project.org/doc/1/reference/generated/kwant.lattice.Polyatomic#kwant.lattice.Polyatomic.shape

Let me know if the documentation or the error message are unclear then.

Best,
Anton

On Mon, 16 Dec 2019 at 10:09, tavakkolidjawad  wrote:
>
>
> Hello everyone
>
> I want to create a lead for a 3D structure, but when I run the program, I get 
> the following error:
>
> ###No sites close to (0, 0, 14) are inside the desired shape.###
>
> Please help me find the mistake.
>
>
>
> Thanks in advance for your reply.
>
>
> 
>
> import kwant
> from matplotlib import pyplot
>
>
> lat = kwant.lattice.general([(0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0)],
> [(0, 0, 0), (0.25, 0.25, 0.25)], name=['a1', 'a2'])
> a1, a2 = lat.sublattices
>
>
> def make_cuboid(a=10, b=10, c=15):
> def cuboid_shape(pos):
> x, y, z = pos
> return 0 <= x < a and 0 <= y < b and 0 <= z < c
>
>
>
> syst = kwant.Builder()
> syst[a1.shape(cuboid_shape, (0, 0, 0))] = 1 ## the same on-site foar a1 and b1
> syst[a2.shape(cuboid_shape, (0, 0, 0))] = 1
> syst[lat.neighbors()] = 1
>
>
> def lead_shape(pos):
> x,y,z = pos
> return 0 <= x < a and 0 <= y < b and 14 <= z < c
>
> lead = kwant.Builder(kwant.TranslationalSymmetry([1,0,0]))
> lead[lat.shape(lead_shape, (0, 0, 14))] = 1
> lead[lat.neighbors()] = 1
>
>
> return syst,lead
>
> def main():
>
> syst,lead = make_cuboid()
>
> kwant.plot(syst)
>
> syst = syst.finalized()
>
> syst = make_cuboid(a=1.1, b=1.1, c=1.1)
>
> def family_colors(site):
> return 'r' if site.family == a1 else 'g'
>
> syst.attach_lead(lead)
>
>
> kwant.plot(syst, site_size=0.18, site_lw=0.01, hop_lw=0.05,
> site_color=family_colors)
>
>
>
> if __name__ == '__main__':
> main()
>
> #

Re: [Kwant] Bound States in a JJ

[Anton Akhmerov]

To add to the discussion, Joe Weston recently implemented a cleaned up
version of the bound state solver, available here:
https://gitlab.kwant-project.org/kwant/kwant/merge_requests/320

If you don't want to install the development version of Kwant, just
grab the boundstate.py file from that merge request.

Best,
Anton

On Tue, 10 Dec 2019 at 09:38, Xavier Waintal  wrote:
>
> Dear Denise,
>
> Slowly decaying bound states can happen whenever the bound state energy is 
> close to the opening
> of a channel (e.g. close to the gap in your case). To move it away from the 
> gap, you can put a phase difference across the superconductors.
>
> As a sanity check you can take a finite system and perform an exact 
> diagonalisation (for several sizes that include a finite portion of the 
> leads). This is imprecise when the bound states are slowly decaying, but it 
> should allow you to spot obvious problems.
>
> Note that the algorithm is still experimental. As it requires some root 
> finding for a function f(E)=0, it is always a good idea to plot the function 
> f(E) to see if the algorithm has missed a root.
>
> Hope this helps.
>
> Best regards,
> Xavier
>
>
> Le 9 déc. 2019 à 19:27, Denise Puglia  a écrit :
>
> Dear All,
>
> Has anyone tried to calculate bound states in a JJ? I used the package by 
> Benoit Rossignol available at: 
> https://gitlab.kwant-project.org/kwant/boundstate .
> I am sending attached the code I used to simulate it but I am having some 
> difficulties interpreting the results. First, the algorithm could not find 
> bound states, i.e., psi_tot returned by the solver is an empty array, for 
> N<5. Second, it returns a total wavefunction for larger N, however the 
> wavefunction oscillates between each site and decays very slowly in the 
> leads. In the original proposal of this method, Istas et al 
> (arXiv:1711.08250) solves the wavefunction of quantum billiard (Fig 5) in 
> which this does not seem to happen. Has anyone found something similar? I do 
> know if it a mistake in the code, an attribute of discretization or something 
> else.
>
> Scanning over N=[1,2,3,5,10,20] and ,u=[2, Delta, 0, -2] the algorithm 
> returns the following non-zero wavefunctions:
> SNS junction with N=10, S=20 and mu=2
> 
> SNS junction with N=20, S=20 and mu=2
> 
> SNS junction with N=5, S=20 and mu=-2 (in this case the amplitude of the 
> bound states decays to zero in the lead)
> 
>
> SNSNS junction with 20/10/10/10/20 sites and mu=2. The wavefunction in the N 
> region is similar to the SNS junction but there seems to be no decay in the 
> middle S region.
> 
> and of course it's mirrored version:
> 
>
> I appreciate any comments on this subject.
>
> Best regards,
> Denise
>
> <_common.py>
>
>

Re: [Kwant] Convert plot of system to .tex file

[Anton Akhmerov]

Hi Alexandre,

I think your best bet is to use tikzplotlib with Kwant plots. Kwant
plotting functions create pyplot figures, so that should work.

Best,
Anton

On Sat, 19 Oct 2019 at 18:59,  wrote:
>
> Hello,
>
>
>
> My name is Alexandre and I would like to know if the file produced by the line
>
> figLattice = kwant.plot( sys, file = 'test.pdf' )
>
> could give the corresponding plot in .tex file (the same way as tikzplotlib 
> works) so that it can be plotted using tikzpicture for example ?
>
> If not, is there a way to get rid of all the "useless" white space around the 
> figure ?
>
>
>
> Thank you and have a good day,
> Alexandre

Re: [Kwant] Plotting energy as a function of magnetic field in 3D.

[Anton Akhmerov]

Dear Naveen,

Yes, the code supports 3D systems, see the last section of the tutorial.

Best,
Anton

On Fri, 13 Sep 2019 at 10:00, Naveen Yadav  wrote:
>
> dear sir,
>
> I have tried it for 3D BHZ model but it is not working. Does it work for only 
> 2D system.
>
>
> On Thu, Sep 12, 2019 at 1:54 PM Anton Akhmerov  
> wrote:
>>
>> Great,
>>
>> Just to add a concluding remark: Landau fan requires a continuum
>> approximation of the Hamiltonian. If you start with a tight-binding
>> Hamiltonian you'd get fractal spectrum and the Hofstadter butterfly
>> instead.
>>
>> Happy Kwanting,
>> Anton
>>
>> On Thu, 12 Sep 2019 at 10:09, Naveen Yadav  wrote:
>> >
>> > The second problem is solved. I got the landau fan for BHZ model.
>> > Thank you for the support.
>> >
>> >
>> > On Thu, Sep 12, 2019, 12:46 Naveen Yadav  wrote:
>> >>
>> >> Dear sir,
>> >>
>> >> I have updated KWANT, but it shows the AttributeError: module 
>> >> 'kwant.continuum' has no attribute 'discretize_landau'. And in browser 
>> >> also, it is showing the same error. I have downloaded the 
>> >> landau_levels.py file and put it into the continuum folder. But it is not 
>> >> working.
>> >>
>> >>
>> >> On Wed, Sep 11, 2019, 23:47 Naveen Yadav  wrote:
>> >>>
>> >>> Dear sir,
>> >>>
>> >>> That is exactly what I am looking for.
>> >>> But in my case Hamiltonian is not polynomial in k. It contains Sine and 
>> >>> Cosine terms. It's a tight binding Hamiltonian having coupling  terms 
>> >>> like sigma_x *Sin(kx)+ sigma_y *sin(ky) and mass term in the 
>> >>> trignometric form. So, can I proceed by writing the trignometric terms 
>> >>> in some lower order polynomial terms? Does that make sense?
>> >>> Please make some comment regarding this.
>> >>> Thank you very much.
>> >>>
>> >>>
>> >>>
>> >>>
>> >>>
>> >>>
>> >>>
>> >>>
>> >>>
>> >>> Naveen
>> >>> Department of Physics & Astrophysics
>> >>> University of Delhi
>> >>> New Delhi-110007
>> >>>
>> >>> On Wed, Sep 11, 2019, 22:18 Anton Akhmerov  
>> >>> wrote:
>> >>>>
>> >>>> Dear Naveen,
>> >>>>
>> >>>> If you are dealing with a continuum Hamiltonian (so a polynomial in
>> >>>> k-space), then there is a recent addition to Kwant, that allows to
>> >>>> compute Landau levels. Please check out if this tutorial is what you
>> >>>> are looking for:
>> >>>> https://kwant-project.org/doc/dev/tutorial/magnetic_field#adding-magnetic-field
>> >>>> (if you click the "activate thebelab" button, you can also play around
>> >>>> with the code in your browser).
>> >>>>
>> >>>> If that suits your needs, you'd need to either install a development
>> >>>> version of Kwant or just get this file:
>> >>>> https://gitlab.kwant-project.org/kwant/kwant/blob/master/kwant/continuum/landau_levels.py
>> >>>>
>> >>>> Let me know if that answers your question,
>> >>>> Anton
>> >>>>
>> >>>> On Wed, 11 Sep 2019 at 18:39, Naveen Yadav  
>> >>>> wrote:
>> >>>> >
>> >>>> > Dear sir,
>> >>>> >
>> >>>> > I understood that this code is off no use. The leads are useless here.
>> >>>> > Actually, I want to plot the Landau fan. Can KWANT  do the job here?
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> >
>> >>>> > Naveen
>> >>>> > Department of Physics & Astrophysics
>> >>>> > University of Delhi
>> >>>> > New Delhi-110007
>> >>>> >
>> >>>> > On Mon, Sep 9, 2019, 00:50 Abbout Adel  wrote:
>> >>>> >>
>> >>>> >> Dear Naveen,
>>

Re: [Kwant] Plotting energy as a function of magnetic field in 3D.

[Anton Akhmerov]

Dear Naveen,

If you are dealing with a continuum Hamiltonian (so a polynomial in
k-space), then there is a recent addition to Kwant, that allows to
compute Landau levels. Please check out if this tutorial is what you
are looking for:
https://kwant-project.org/doc/dev/tutorial/magnetic_field#adding-magnetic-field
(if you click the "activate thebelab" button, you can also play around
with the code in your browser).

If that suits your needs, you'd need to either install a development
version of Kwant or just get this file:
https://gitlab.kwant-project.org/kwant/kwant/blob/master/kwant/continuum/landau_levels.py

Let me know if that answers your question,
Anton

On Wed, 11 Sep 2019 at 18:39, Naveen Yadav  wrote:
>
> Dear sir,
>
> I understood that this code is off no use. The leads are useless here.
> Actually, I want to plot the Landau fan. Can KWANT  do the job here?
>
>
>
>
>
>
>
>
>
>
>
> Naveen
> Department of Physics & Astrophysics
> University of Delhi
> New Delhi-110007
>
> On Mon, Sep 9, 2019, 00:50 Abbout Adel  wrote:
>>
>> Dear Naveen,
>>
>> If your concern is the program which is slow, that is not an issue since it 
>> takes just few minutes.
>> Now, if you are talking about the result, I want to be sure that you notice 
>> that your system is not infinite as you claim in your email.
>> You can check that by adding extra cells from the lead" 
>> syst.attach_lead(lead, add_cells=10)
>> Actually, in your case, the presence of the leads is useless since at the 
>> end, you are just diagonalizing the Hamiltonian of the central system.
>> If you want to study an infinite system in x and y, you need to look at the 
>> module "wraparound" and the example of graphene that is in the archive of 
>> kwant.
>> For the magnetic field, you can use the Pierls substitution. check for 
>> example this paper [1]
>>
>> You can also think about the use of continuous Hamiltonian in kwant. You may 
>> find it very useful [2]
>> I hope this helps.
>>
>> Regards,
>> Adel
>>
>>
>> [1]  https://arxiv.org/pdf/1601.06507.pdf
>> [2] https://kwant-project.org/doc/1/tutorial/discretize
>>
>> On Sun, Sep 8, 2019 at 6:16 PM Naveen Yadav  wrote:
>>>
>>> Dear Sir,
>>> Thanks for the tips. As you told, I have tried in other way also but I am 
>>> getting the same result which are very tedious. I don't know where is fault.
>>> Now the code looks like
>>>
>>> import kwant
>>> import scipy.sparse.linalg as sla
>>> import matplotlib.pyplot as plt
>>> import tinyarray
>>> import numpy as np
>>> from numpy import cos, sin, pi
>>> import cmath
>>> from cmath import exp
>>>
>>> sigma_0 = tinyarray.array([[1, 0], [0, 1]])
>>> sigma_x = tinyarray.array([[0, 1], [1, 0]])
>>> sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
>>> sigma_z = tinyarray.array([[1, 0], [0, -1]])
>>>
>>>
>>> def make_system(a=1, L=30, W=10, H=10, t=1.0, t_x=1.0, t_y=1.0, t_z=1.0, 
>>> lamda=0.1, beta=1.05):
>>> def onsite(site):
>>> return (t_z * cos(beta) + 2 * t) * sigma_z
>>>
>>> def hoppingx(site0, site1):
>>> return (-0.5 * t * sigma_z - 0.5 * 1j * t_x * sigma_x)
>>>
>>> def hoppingy(site0, site1):
>>> return -0.5 * t * sigma_z - 0.5 * 1j * t_y * sigma_y
>>>
>>> def hoppingz(site0, site1, B):
>>> y = site1.pos[1]
>>> return (-0.5 * t_z * sigma_z - 0.5 * 1j * lamda * sigma_0) * exp(2 
>>> * pi * 1j * B * a * (y-40))
>>>
>>>
>>> syst = kwant.Builder()
>>> lat = kwant.lattice.cubic(a)
>>> syst[(lat(z, y, x) for z in range(H) for y in range(W) for x in 
>>> range(L))] = onsite
>>> syst[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz
>>> syst[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy
>>> syst[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx
>>>
>>> lead1=kwant.Builder(kwant.TranslationalSymmetry((0,-a,0)))
>>> lead1[(lat(z,y,x)  for z in range(H)for y in range(W)for x in 
>>> range(L))]=onsite
>>> lead1[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz
>>> lead1[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy
>>> lead1[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx
>>>
>>> syst.attach_lead(lead1)
>>> syst.attach_lead(lead1.reversed())
>>>
>>> lead2=kwant.Builder(kwant.TranslationalSymmetry((-a,0,0)))
>>> lead2[(lat(z,y,x)  for z in range(H)for y in range(W)for x in 
>>> range(L))]=onsite
>>> lead2[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz
>>> lead2[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy
>>> lead2[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx
>>>
>>> syst.attach_lead(lead2)
>>> syst.attach_lead(lead2.reversed())
>>> syst = syst.finalized()
>>> return syst
>>>
>>> def analyze_system(syst, Bfields):
>>> syst = make_system()
>>> kwant.plot(syst)
>>> energies = []
>>> for B in Bfields:
>>> #print(B)
>>> ham_mat = syst.hamiltonian_submatrix(params=dict(B=B), sparse=True)
>>> 

Re: [Kwant] Wavefunction for 2D NISIN

[Anton Akhmerov]

Hi Hanifa,

This is certainly possible. It's up to you for which wave function you
compute the density. You can compute the eigenstates and compute the
expectation all the same.

Best,
Anton

On Mon, 26 Aug 2019 at 18:46, Tidjani, Hanifa  wrote:
>
> Dear Anton,
>
>
>
> Thank you for your reply.
>
>
>
> This does indeed work, however I am concerned that I am driving it at a 
> certain energy and then plotting the wave function. Is there a way to plot 
> the wave function of a particular eigenstate rather than at a particular 
> energy as is done in tutorial 2.4?
>
>
>
> Kind regards,
>
>
>
> Hanifa
>
>
>
> 
> From: Anton Akhmerov 
> Sent: Monday, August 26, 2019 9:16:59 AM
> To: Tidjani, Hanifa 
> Cc: kwant-discuss@kwant-project.org 
> Subject: Re: [Kwant] Wavefunction for 2D NISIN
>
> Hi Hanifa,
>
> The mismatch between the eigenvector size and the number of sites in
> the system is because your system has more than one degree of freedom
> per site. Please see how to deal with such cases in the tutorial:
> https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fkwant-project.org%2Fdoc%2F1%2Ftutorial%2Foperatorsdata=01%7C01%7Chanifa.tidjani%40kcl.ac.uk%7C4d22e15c6b5845a7382908d729fdcb8f%7C8370cf1416f34c16b83c724071654356%7C0sdata=bLzi8aLHdVqCfwSycyYE4ZLXhpsJeZR15iUHFF0pZgE%3Dreserved=0
>
> Best,
> Anton
>
> On Mon, 26 Aug 2019 at 10:08, Tidjani, Hanifa  
> wrote:
> >
> >
> >
> > I am trying to plot the wavefunction for a 2D NISIN Hamiltonian however I 
> > am unable to get a result, this is my code:
> >
> >
> >
> >
> >
> > import kwant
> >
> > import tinyarray as ta
> >
> > import numpy as np
> >
> > import matplotlib.pyplot as plt
> >
> > from tqdm import tqdm
> >
> > import winsound
> >
> > import os
> >
> > from datetime import datetime
> >
> > import scipy.sparse.linalg as sla
> >
> >
> >
> >
> >
> > sig_0 = np.identity(4)
> >
> > s_0 = np.identity(2)
> >
> > s_z = np.array([[1., 0.], [0., -1.]])
> >
> > s_x = np.array([[0., 1.], [1., 0.]])
> >
> > s_y = np.array([[0., -1j], [1j, 0.]])
> >
> >
> >
> > tau_z = ta.array(np.kron(s_z, s_0))
> >
> > tau_x = ta.array(np.kron(s_x, s_0))
> >
> > tau_y = ta.array(np.kron(s_y, s_0))
> >
> > sig_z = ta.array(np.kron(s_0, s_z))
> >
> > tau_zsig_x = ta.array(np.kron(s_z, s_x))
> >
> > tau_zsig_y = ta.array(np.kron(s_z, s_y))
> >
> >
> >
> > # Lorentzian definition
> >
> >
> >
> >
> >
> > def lorentzianxy(x, y, ind,y0):
> >
> > gam = 3
> >
> > L = params["L"]
> >
> > #if 0 >
> > return gam ** 2 /( gam ** 2 + (x - ind) ** 2 + (y - y0) ** 2)
> >
> ># else :
> >
> > #return 0
> >
> >
> >
> > # Define potential
> >
> >
> >
> > def potential(site, pot,ind,y0):
> >
> > (x, y) = site.pos
> >
> > return pot * lorentzianxy(x,y,ind,y0)
> >
> >
> >
> >
> >
> > # Make system
> >
> >
> >
> > def makeNISIN2D(params):
> >
> >
> >
> > # List all the params
> >
> > W=params["W"]
> >
> > L=params["L"]
> >
> > pot = params["pot"]
> >
> > ind = params["ind"]
> >
> > mu=params["mu"]
> >
> > B=params["B"]
> >
> > Delta=params["Delta"]
> >
> > alpha=params["alpha"]
> >
> > t=params["t"]
> >
> > barrier = params["barrier"]
> >
> > Smu = params["Smu"]
> >
> >
> >
> >
> >
> >
> >
> > def lorentzianxy(x, y, ind,y0):
> >
> > gam = 3
> >
> > L = params["L"]
> >
> > #if 0 >
> > return gam ** 2 /( gam ** 2 + (x - ind) ** 2 + (y - y0) ** 2)
> >
> ># else :
> >
> > #return 0
> >
> >
> >
> > # Define potential
> >
> >
> >
> > def potential(site, pot,ind,y0):
> >
> > (x, y) = site.pos
> >
> > return pot * lorentzianxy(x,y,ind,y0)
> >
> > #if y < lorentzian(x,ind)*4:
> >
>

Re: [Kwant] Wavefunction for 2D NISIN

[Anton Akhmerov]

Hi Hanifa,

The mismatch between the eigenvector size and the number of sites in
the system is because your system has more than one degree of freedom
per site. Please see how to deal with such cases in the tutorial:
https://kwant-project.org/doc/1/tutorial/operators

Best,
Anton

On Mon, 26 Aug 2019 at 10:08, Tidjani, Hanifa  wrote:
>
>
>
> I am trying to plot the wavefunction for a 2D NISIN Hamiltonian however I am 
> unable to get a result, this is my code:
>
>
>
>
>
> import kwant
>
> import tinyarray as ta
>
> import numpy as np
>
> import matplotlib.pyplot as plt
>
> from tqdm import tqdm
>
> import winsound
>
> import os
>
> from datetime import datetime
>
> import scipy.sparse.linalg as sla
>
>
>
>
>
> sig_0 = np.identity(4)
>
> s_0 = np.identity(2)
>
> s_z = np.array([[1., 0.], [0., -1.]])
>
> s_x = np.array([[0., 1.], [1., 0.]])
>
> s_y = np.array([[0., -1j], [1j, 0.]])
>
>
>
> tau_z = ta.array(np.kron(s_z, s_0))
>
> tau_x = ta.array(np.kron(s_x, s_0))
>
> tau_y = ta.array(np.kron(s_y, s_0))
>
> sig_z = ta.array(np.kron(s_0, s_z))
>
> tau_zsig_x = ta.array(np.kron(s_z, s_x))
>
> tau_zsig_y = ta.array(np.kron(s_z, s_y))
>
>
>
> # Lorentzian definition
>
>
>
>
>
> def lorentzianxy(x, y, ind,y0):
>
> gam = 3
>
> L = params["L"]
>
> #if 0
> return gam ** 2 /( gam ** 2 + (x - ind) ** 2 + (y - y0) ** 2)
>
># else :
>
> #return 0
>
>
>
> # Define potential
>
>
>
> def potential(site, pot,ind,y0):
>
> (x, y) = site.pos
>
> return pot * lorentzianxy(x,y,ind,y0)
>
>
>
>
>
> # Make system
>
>
>
> def makeNISIN2D(params):
>
>
>
> # List all the params
>
> W=params["W"]
>
> L=params["L"]
>
> pot = params["pot"]
>
> ind = params["ind"]
>
> mu=params["mu"]
>
> B=params["B"]
>
> Delta=params["Delta"]
>
> alpha=params["alpha"]
>
> t=params["t"]
>
> barrier = params["barrier"]
>
> Smu = params["Smu"]
>
>
>
>
>
>
>
> def lorentzianxy(x, y, ind,y0):
>
> gam = 3
>
> L = params["L"]
>
> #if 0
> return gam ** 2 /( gam ** 2 + (x - ind) ** 2 + (y - y0) ** 2)
>
># else :
>
> #return 0
>
>
>
> # Define potential
>
>
>
> def potential(site, pot,ind,y0):
>
> (x, y) = site.pos
>
> return pot * lorentzianxy(x,y,ind,y0)
>
> #if y < lorentzian(x,ind)*4:
>
> #
>
> #return (tru_pot)
>
> #else:
>
> #return 0
>
> #
>
> def Straightpotential(site, pot,ind):
>
> (x, y) = site.pos
>
> if x == ind:
>
> return pot
>
> else:
>
> return 0
>
>
>
> def onsiteSc(site, pot,ind,y0):
>
> #(x,y)=site.pos
>
> #L=params["L"]
>
># if x>L/2:
>
>return (4 * t - Smu) * tau_z + B * sig_z + Delta * tau_x - 
> potential(site, pot,ind,y0)*tau_z
>
># else:
>
>#return (4 * t ) * tau_z + B * sig_z + Delta * tau_x - potential(site, 
> pot,ind,y0)*tau_z
>
>
>
> def onsiteNormal(site, pot,ind,y0):
>
> return (4 * t - mu) * tau_z #- potential(site, pot,ind,y0)*tau_z
>
> def onsiteBarrier(site, pot,ind,y0):
>
> return (4 * t - mu + barrier) * tau_z #- potential(site, 
> pot,ind,y0)*tau_z
>
>
>
> def hop_x(site, pot, ind,y0):
>
>
>
> return -t * tau_z + 0.5j * alpha * tau_zsig_y
>
>
>
> def hop_y(site, pot, ind,y0):
>
>
>
>return -t * tau_z - .5j * alpha * tau_zsig_x
>
>
>
>
>
>
>
> # On each site, electron and hole orbitals.
>
> lat = kwant.lattice.square(norbs=4)
>
> syst = kwant.Builder()
>
>
>
> # S
>
> syst[(lat(i, j) for i in range(1,L-1) for j in range(W))] = onsiteSc
>
>
>
> barrierLen = 1;
>
> # I's
>
> syst[(lat(i, j) for i in range(barrierLen) for j in range(W))] = 
> onsiteBarrier
>
> syst[(lat(i, j) for i in range(L-barrierLen, L)for j in range(W))] = 
> onsiteBarrier
>
> syst[kwant.builder.HoppingKind((1, 0), lat, lat)]= hop_x
>
> syst[kwant.builder.HoppingKind((0, 1), lat, lat)]= hop_y
>
>
>
> # N's
>
> lead = kwant.Builder(kwant.TranslationalSymmetry((-1,0)),
>
>  conservation_law=-tau_z#, particle_hole=tau_y
>
>  )
>
> lead[(lat(0, j) for j in range(W))] = onsiteNormal
>
> lead[kwant.builder.HoppingKind((1, 0), lat, lat)]= hop_x
>
> lead[kwant.builder.HoppingKind((0, 1), lat, lat)]= hop_y
>
>
>
> syst.attach_lead(lead)
>
> syst.attach_lead(lead.reversed())
>
>
>
> return syst
>
>
>
>
>
> def sorted_eigs(ev):
>
> evals, evecs = ev
>
> evals, evecs = map(np.array, zip(*sorted(zip(evals, evecs.transpose()
>
> return evals, evecs.transpose()
>
>
>
> params = dict(mu=0.3, Delta=.1, alpha=.8, t=1.0, barrier = 2.0, pot = 0.0,W = 
> 5, L = 30, ind = 0, B = 0.3, Smu=0.00,y0=0)
>
>
>
> sys = makeNISIN2D(params)
>
> sys = sys.finalized()
>
>
>
> # Calculate the wave functions in the system.
>
> ham_mat = 

Re: [Kwant] Hopping settings

[Anton Akhmerov]

Dear Hosein,

You are creating a system with essentially all possible hoppings.
These will take a considerable computational cost regardless of what
you do. Try limiting the number of hoppings to only a few nearest
neighbors per site.

Best,
Anton

On Wed, 14 Aug 2019 at 09:40, Khani Hosein  wrote:
>
> Dear developers,
> I constructed Lat_A and Lat_B for a system.
> A1,A2=Lat_A.sublattices
> B1,B2=Lat_B.sublattices
> I want to set hoppings between Lat_A and Lat_B, and the hoppings are 
> dependent on the site position.Every atom of Lat_A is coupled to all the 
> atoms in Lat_B.
>
> Now I have tried it like this :
> def coupling(site1,site2):
> xp1,yp1=site1.pos
> xp2,yp2=siye2.pos
> return t*sqrt((xp1-xp2)**2+(yp1-yp2)**2)/d
>
> for i in range(100):
> for j in range(100):
> sys[kwant.builder.HoppingKind(( i, j), A1, B1)] = coupling
> sys[kwant.builder.HoppingKind(( i, j), A1, B2)] = coupling
> sys[kwant.builder.HoppingKind(( i, j), A2, B1)] = coupling
> sys[kwant.builder.HoppingKind(( i, j), A2, B2)] = coupling
>
> It can not work because it take too long to do the calculations. I am 
> considering a bilayer system but rotations are considered, for example a 
> twisted bilayer graphene. Could you please help me to solve this problem?
> Regards,
> Hosein Khani

Re: [Kwant] How to feed a portion of a finite system to kwant.plotter.interpolate_density

[Anton Akhmerov]

Hi Hadi,

The intended way would be to extend the density vector with np.nan, so
that it has a correct length. However since you mention performance,
I'm not sure if this would achieve what you are looking for.
Inspecting the source of interpolate_density [1], we see that the only
way the system is used is to get

sites = np.array([s.pos for s in syst.sites])

Therefore we can mock the necessary properties of the system by using
the following code:

from types import SimpleNamespace
# make a fake system that only contains the correct sites and isn't
actually a system
syst_slice = SimpleNamespace(sites=(site for site in syst.sites if where(site)))

It is certainly a hack, so a more reliable solution would be to take
over the bits of the source of interpolate_density and the low level
function that it uses [2] and re-purpose for your needs.

[1]: 
https://gitlab.kwant-project.org/kwant/kwant/blob/v1.4.1/kwant/plotter.py#L1793-1872
[2]: 
https://gitlab.kwant-project.org/kwant/kwant/blob/v1.4.1/kwant/plotter.py#L1637-1701

On Fri, 16 Aug 2019 at 10:03, Hadi Zahir  wrote:
>
> Dear Kwant developers,
>
> I have a finalized closed system ''sys_f'' for which I computed the local 
> charge densities "ch_densities", which is basically a vector of 
> length=len(sys_f.sites).
>
> For complexity reason,  I'd like to feed a slice of sys_f  to the 
> interpolate_density() as:
>
> >> kwant.plotter.interpolate_density(syst=sys_p, density=ch_densities_p)
>
> where sys_p is a slice of sys_f and ch_densities_p is the corresponding 
> density vector.
> Could you please let me know how to do it?
>
> Best regards,
> Hadi

Re: [Kwant] 3D annulus cylinder

[Anton Akhmerov]

Dear Naveen,

These inwards- and outwards- facing leads wouldn't have translation
invariance. This makes one unable to compute a mode decomposition in
such a geometry, and therefore computing the conductance becomes
extremely difficult. Right now Kwant only implements the algorithms
that assume translationally invariant leads.

Best,
Anton

On Tue, 13 Aug 2019 at 14:15, Naveen Yadav  wrote:
>
> Dear sir,
>
> syst.attach_lead(lead, origin=lat(0, 0, 0))  # lat(0, 0, 0) is in the hole of 
> the annulus
> This is okay. But I want to create leads in the radial direction, suppose X 
> is the width of cylinder,  Y is circumference and Z is the Difference in 
> outer and inner radii. So, I want to create leads wrapped around Y, for inner 
> circumference lead should directed towards origin and for outer circle 
> directed away from the origin. Code for creating annulus geometry is given 
> below-
>
> import kwant
> import scipy.sparse.linalg as sla
> import matplotlib.pyplot as plt
> import tinyarray
> import numpy as np
> from numpy import cos, sin, pi
> import cmath
> from cmath import exp
>
> sigma_0 = tinyarray.array([[1, 0], [0, 1]])
> sigma_x = tinyarray.array([[0, 1], [1, 0]])
> sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
> sigma_z = tinyarray.array([[1, 0], [0, -1]])
>
>
> def make_system(a=1, L=22, r_in=22, r_out=30, t=1.0, t_x=1.0, t_y=1.0, 
> t_z=1.0, lamda=0.2, beta=1.05, phi_uc = 0.0078):
> # ring shape
> def ring(pos):
> (z, y, x) = pos
> rsq = y ** 2 + z ** 2
> return r_in ** 2 <= rsq <= r_out ** 2 and x in range (L)
>
> def onsite(site):
> return (t_z * cos(beta) + 2 * t) * sigma_z
>
> def hoppingx(site0, site1):
> return (-0.5 * t * sigma_z - 0.5 * 1j * t_x * sigma_x)
>
> def hoppingy(site0, site1):
> return -0.5 * t * sigma_z - 0.5 * 1j * t_y * sigma_y
>
> def hoppingz(site0, site1):
> y = site1.pos[1]
> return (-0.5 * t_z * sigma_z - 0.5 * 1j * lamda * sigma_0) * exp(2 * 
> pi * 1j * phi_uc * a * (y-40))
>
>
> syst = kwant.Builder()
> lat = kwant.lattice.cubic(a, norbs=2)
> syst[lat.shape(ring, (0, r_in+1, 0))] = onsite
> syst[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz
> syst[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy
> syst[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx
>
> lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0, 0)))
>
> lead[lat.shape(ring, (0, r_in+1, 0))] = onsite
> lead[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz
> lead[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy
> lead[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx
>
> syst.attach_lead(lead)
> syst.attach_lead(lead, origin=lat(0,0,0))
> return syst
>
> def analyze_system():
> syst = make_system()
> fig = plt.figure()
> ax = kwant.plot(syst)
> ax.savefig('sys2.png',dpi=200)
>
> def main():
> syst = make_system()
> analyze_system()
> main()
>
> On Tue, Aug 13, 2019 at 3:18 PM Joseph Weston  
> wrote:
>>
>> Hi,
>>
>> Dear sir,
>>
>> Could we attach circular leads to the inner and outer circle of annulus 
>> geometry in 3D? Please suggest me if there is a way to do that.
>>
>> What do you mean by circular leads? Do you mean leads with a circle 
>> cross-section (i.e. a semi-infinite cylinder lead)? If so then all you need 
>> to do is create a lead with a circular cross-section uses 'lat.shape' in a 
>> similar way to how you created the scattering region. You can attach leads 
>> to the interior of the annulus by specifying the parameter 'origin' to be a 
>> site in the interior of the annulus when calling 'attach_lead' (see the 
>> documentation [1]). e.g.:
>>
>> syst.attach_lead(lead, origin=lat(0, 0))  # lat(0, 0) is in the hole of 
>> the annulus
>>
>>
>> Happy Kwanting,
>>
>>
>> Joe
>>
>>
>> [1]: 
>> https://kwant-project.org/doc/1/reference/generated/kwant.builder.Builder#kwant.builder.Builder.attach_lead
>
>
>
> --
>
>
> With Best Regards
> NAVEEN YADAV
> Ph.D Research Scholar
> Deptt. Of Physics & Astrophysics
> University Of Delhi.

Re: [Kwant] Group velocites in the leads

[Anton Akhmerov]

Dear Luca,

If done correctly, both approaches lead to the same answer, but Kwant
utilizes the second one.

Best regards,
Anton Akhmerov

On Mon, 29 Jul 2019 at 16:57, Luca Lepori  wrote:
>
> Dear Kwant developers,
>
> I write concerning the group velocites
> of the modes in the leads.
> In particular, I wonder how Kwant calculates them.
>
> I see two possibilities (leads along x):
> 1) knowing the dispersion along x, from the usual definition v(kx) =
> \partial E(kx) / \partial kx;
> 2) working equivalently in the real space basis, and in particular on
> a unitary cell
> of a lead, one can calculate vx also as:
> vx = <\psi |i [H, \hat{x}] | \psi>,
> being \psi the wavefunction of the mode and
> H the Hamiltonian respectively, both defined on the unitary cell of
> the lead, and \hat{x}
> the position operator.
> The second approach can be useful also to estimate the expectation
> value of the velocities in the direction orthogonal to the leads (say
> when close boundary conditions are assumed).
>
> Does Kwant utilizes directly one of these two strategies ?
>
> Thank you very much and best regards
>
> L. L.

Re: [Kwant] Periods for kwant.TranslationalSymmetry

[Anton Akhmerov]

Hi Barış, everyone,

Barıș, you aim to have a lead with a different (rotated) lattice,
compared to the scattering region. This is possible, but requires
extra manual work (see
https://kwant-project.org/doc/1/tutorial/faq#how-to-use-different-lattices-for-the-scattering-region-and-a-lead).

I am, however, not sure if this is really what you want. In most
mesoscopic systems, such a lattice change doesn't occur naturally, and
if it does, then it does not have a very profound effect.

Best,
Anton

On Mon, 15 Jul 2019 at 21:10, Barış Özgüler  wrote:
>
> Hi Ousmane et al,
>
> I actually want the block of unit cells to be perpendicular to the scattering 
> region as in the Wanted_System.png, attached. Note that it is a 
> representative plot, I didn't produce it via Kwant but I want to attach the 
> red lead to the black region in that direction, not like in 
> Unwanted_System.png, which is also attached and produced via the code:
>
> ##
>
> import kwant
> import numpy as np
> import matplotlib.pyplot as plt
> from types import SimpleNamespace
>
> def make_system(a, W, L):
>
> def shape(pos):
>   (x, y) = pos
>
>   return (0 <= y <= W/2 and -x/2 <= y <= (-x + L)/2) # or (-W <= y <= 
> 0 and 0 <= x <= L)
>
> def onsite(site, par):
> return 4 * par.t - par.mu
>
> def hopx(par):
> return -par.t
>
> def hopy(par):
> return -par.t
>
>
> # lead
> lat = kwant.lattice.square(a, norbs=1)
> syst = kwant.Builder()
> syst[lat.shape(shape, (0, 0))] = onsite
> syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
> syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy
> lead = kwant.Builder(kwant.TranslationalSymmetry((-a, -2*a)))
>
> def lead_shape(pos):
> (x, y) = pos
> return -W <= x <= 0
>
> def lead_hopx(site1, site2, par):
> return -par.t
>
> def lead_onsite(site, par):
> return 4 * par.t - par.mu
>
> lead[lat.shape(lead_shape, (-2*a, -a))] = lead_onsite
> lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = lead_hopx
> lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy
> syst.attach_lead(lead)
>
>
>
> syst = syst.finalized()
> return syst
>
>
> syst = make_system(a=5, W=100, L=100)
> kwant.plot(syst)
>
> ##
>
> I couldn't make the lead rotated in the perpendicular direction to the 
> surface. Is there a way to do it?
>
> Thanks,
>
> Barış

Re: [Kwant] Phlead.py

[Anton Akhmerov]

Hi Naveen,

Please take a look at the new superconductivity tutorial:
https://kwant-project.org/doc/1/tutorial/superconductors

That's the most up to date and systematic way.

Best,
Anton

On Mon, 15 Jul 2019 at 13:14, Naveen Yadav  wrote:
>
> Dear all,
>
> Could anyone provide me the newer version of code phlead.py to calculate the 
> conductance through majorana. The older vesion of code is not working well 
> with latest version of kwant and it shows a lot of errors. Please help me in 
> this regard.
>
> Thank you.
> Best regards
> Naveen
>
>
>
>
>
>
>
>
>
>
>
>
>
> Naveen
> Department of Physics & Astrophysics
> University of Delhi
> New Delhi-110007

Re: [Kwant] bounds on KPM

[Anton Akhmerov]

Hi Antonio,

KPM is only stable when the bounds include the complete Hamiltonian.
Otherwise the decomposition terms diverge. This should still work fine with
a small number of moments, but I don't recommend to try KPM in this way.
I'm actually surprised that it worked with LDOS; could it be that you only
considered a small amount of moments?

Cheers,
Anton

On Mon, Jun 17, 2019 at 1:48 PM Abbout Adel  wrote:

> Dear Antonio,
>
> I notice that you result is infinite:  on the y-axis  you have 1e206 !
> I guess that you are dealing with a singularity. I see that you are using
> the bound [0,0.4]. In graphene, the energy 0 should be avoided in  general
> since some singularities appear there (for example the S matrix is singular
> at E=0). Could you try with [.1,.4] for example?
>
> Try to give a reproducible example (if not confidential) in order to
> check  the problem.
>
> I hope this helps
> Adel
>
>
> On Mon, Jun 17, 2019 at 2:35 PM Antonio Lucas Rigotti Manesco <
> antonio...@usp.br> wrote:
>
>> Hi Adel,
>>
>> I'm doing some tests with graphene in the presence of a spatial-dependent
>> magnetic field. First, I calculated the DOS:
>> spectrum = kwant.kpm.SpectralDensity(fsyst)
>> and LDOS:
>> kwant_op = kwant.operator.Density(fsyst, where=sm_shape, sum=False)
>> local_dos = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op,
>> mean=False)
>> and both calculations lead to the results I expected.
>>
>> Next, I did the same calculation, with lower magnetic field (to reproduce
>> real-world conditions), and I noticed that the required memory needed was
>> way higher. So, I tried to fix a smaller window for KPM bounds:
>> spectrum = kwant.kpm.SpectralDensity(fsyst, bounds=(0.,0.4))
>> and
>> local_dos = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op,
>> mean=False, bounds=(0.,0.4))
>>
>> The fact is: LDOS was similar in both calculations. On the other hand,
>> the resulting DOS looks like this:
>> [image: image.png]
>> which do not remind me about the Landau levels of the first calculation
>> at all.
>>
>> So, I'm not sure if I can rely on the LDOS result with limited energy
>> window, since DOS is different from the expected. Also, I wonder if it is
>> possible to get the correct DOS with limited energy window.
>>
>> Best,
>> Antonio
>>
>> Em dom, 16 de jun de 2019 às 05:11, Abbout Adel 
>> escreveu:
>>
>>> Dear Antonio,
>>>
>>> Without an example, your problem may be not very clear. Could you
>>> provide an example?
>>>
>>> I think that you are talking about the fact that the values that you get
>>> when you change the  length of the interval  are different which is quite
>>> normal.
>>> Put in mind that, the result is the number of appearance of some
>>> energies and that the density is not normalized.
>>>
>>> Regards,
>>> Adel
>>>
>>> On Thu, Jun 13, 2019 at 7:46 PM Antonio Lucas Rigotti Manesco <
>>> antonio...@usp.br> wrote:
>>>
 Hello all,

 I'm trying to do some calculations using KPM method for graphene in a
 somehow large quantum dot (300ax300a, where a is the lattice constant).
 Since LDOS calculations uses lot of memory and I need good energy
 resolution, I redefined my energy bounds to a smaller window, and the
 behavior is similar to what I expected. However, if I do the same for DOS
 calculations, I get something that looks very strange, compared with the
 previous result using the automatically generated bounds. So, my question
 is: are the LDOS results reliable even though the DOS does not looks right?
 (I'm calculating the LDOS for energies about the middle of the new window,
 in case that matters).

 Thanks
 Antonio

>>>
>>>
>>> --
>>> Abbout Adel
>>>
>>
>>
>> --
>> Antônio Lucas Rigotti Manesco
>> PhD fellow - University of São Paulo, Brazil
>>
>
>
> --
> Abbout Adel
>

Re: [Kwant] translational symmetry for non-integer polyatomic basis

[Anton Akhmerov]

Hi Hang Zang,

Can you please explain what you want to achieve? It seems you are
trying to create sites that belong to a lattice, but do not occupy
lattice positions.

Best,
Anton

On Thu, Jun 6, 2019 at 5:01 PM Zang, Hang  wrote:
>
> Dear Kwant developers,
>
> I want to calculate the band structure for a system with the polyatomic basis,
> but since the periodic value is not an integer, there was some problem in 
> finalizing the system. Could you give me some suggestions to solve this 
> problem?
>
> Enclose please find my kwant script.
>
> Thank you very much for your help.
> Hang Zang
>
> ===
> import kwant
> from matplotlib import pyplot
>
> primitive_vectors = [ (1.0,0.0,0.0) , (0.0,1.0,0.0) , (0.0,0.0,1.0) ]
>
> atomic_sites = [ (1.0,1.0,1.0) , (0.0,0.0,0.0) , (-1.0,-1.0,-1.0) ]
>
> lat = kwant.lattice.Polyatomic(primitive_vectors , atomic_sites)
>
> atoms = lat.sublattices
>
> sym_lead = kwant.TranslationalSymmetry((-1.0,0.0,0.0))
>
> lead = kwant.Builder(sym_lead)
>
> # inside the lead
> lead[atoms[0](0.0,0.0,0.0)] = 1
> lead[atoms[1](0.0,0.0,0.0)] = 1
> lead[atoms[2](0.0,0.0,0.0)] = 1
>
> lead[atoms[0](0.0,0.0,0.0),atoms[1](0.0,0.0,0.0)] = -2
> lead[atoms[1](0.0,0.0,0.0),atoms[2](0.0,0.0,0.0)] = -2
>
> # between adjacent leads
> lead[atoms[0](0.0,0.0,0.0),atoms[0](-2.3,0.0,0.0)] = -4
> lead[atoms[1](0.0,0.0,0.0),atoms[1](-2.3,0.0,0.0)] = -4
> lead[atoms[2](0.0,0.0,0.0),atoms[2](-2.3,0.0,0.0)] = -4
>
> lead = lead.finalized()
>
> kwant.plotter.bands(lead,show=False)
> pyplot.xlabel("momentum [(lattice constant)^-1]")
> pyplot.ylabel("energy [t]")
> pyplot.show()
> =

Re: [Kwant] Changing Hamiltonian after finalized()

[Anton Akhmerov]

Hi Pavlo,

- You may make the Hamiltonian have an arbitrary functional dependence
on a parameter say def h(site, mu, x0): return mu * tanh(site.pos[0] -
x0) This is described in this tutorial:
https://kwant-project.org/doc/1/tutorial/spin_potential_shape#spatially-dependent-values-through-functions
- For even more flexibility you may make an arbitrary function a
parameter of your Hamiltonian, e.g. def h(site, mu): return
mu(site.pos)

Best,
Anton

On Thu, May 23, 2019 at 10:14 PM Pavlo Bulanchuk
 wrote:
>
> Dear Kwant developers,
>
> I am trying to calculate conductance of a system with disorder at different 
> chemical potentials. The geometry is constant, so I want to finalize() only 
> once. The question:how do I change onsite potential after finalization?
>
> There is an example from the manual with conductance through a well of 
> variable depth, but if I do it that way, I would generate new onsite 
> potential every time I change chemical potential, which I would like to avoid.
>
> Thank you,
> Pavlo

Re: [Kwant] Error using Qsymm in kwant

[Anton Akhmerov]

Hi Srilok,

Firstly, kwant.qsymm works on builders, so you shouldn't finalize the
system. Furthermore, you shouldn't apply wraparound to the system
(qsymm requires the original symmetry). After you fix these two
things, qsymm will work on your system.

Best,
Anton

On Tue, Mar 19, 2019 at 12:35 AM Srilok Srinivasan  wrote:
>
> Dear kwant developers,
>
> I am currently testing the new feature of kwant 1.4 - the interface with 
> qsymm library.
>
> kag_lat = kwant.lattice.kagome()
> kagome = kwant.Builder(kwant.TranslationalSymmetry(*kag_lat.prim_vecs))
> kagome[kag_lat.shape((lambda pos: True), (0, 0))] = 0
> kagome[kag_lat.neighbors(1)] = 1
> kwant.plot(kagome,site_color=family_color,site_lw=0.01)
> kagome = wraparound(kagome).finalized()
> #dispersion_2D(kag_lat,kagome)
> import kwant.qsymm
> kwant.qsymm.find_builder_symmetries(kagome)
>
> When I try to find the symmetries of the kagome lattice using the 
> find_builder_symmetries I get the following error. Can you please help me fix 
> this?
>
> ---
> TypeError Traceback (most recent call last)
>  in 
>   7 #dispersion_2D(kag_lat,kagome)
>   8 import kwant.qsymm
> > 9 kwant.qsymm.find_builder_symmetries(kagome)
>
> /anaconda3/lib/python3.6/site-packages/kwant/qsymm.py in 
> find_builder_symmetries(builder, momenta, params, spatial_symmetries, 
> prettify, sparse)
> 449
> 450 ham = builder_to_model(builder, momenta=momenta,
> --> 451real_space=False, params=params)
> 452
> 453 # Use dense or sparse computation depending on Hamiltonian size
>
> /anaconda3/lib/python3.6/site-packages/kwant/qsymm.py in 
> builder_to_model(syst, momenta, real_space, params)
> 153 # translation vectors.
> 154 if dim == 0:
> --> 155 proj = np.empty((0, len(list(syst.sites())[0].pos)))
> 156 elif dim < len(list(syst.sites())[0].pos):
> 157 proj, r = la.qr(np.array(periods).T, mode='economic')
>
> TypeError: 'tuple' object is not callable
>
>
> Thanks
>
> Regards,
>
> Srilok Srinivasan
>
>

Re: [Kwant] Kwant 1.4 version

[Anton Akhmerov]

Dear Ali,

You should install Kwant 1.4.0 instead of the pre-release version.

Best,
Anton

On Mon, Mar 18, 2019 at 12:41 PM Ali Asgharpour
 wrote:
>
> Hello
>
> I think I have updated kwant on my laptop since it gives me the version 
> 1.4.0a1 when I write kwant.version.version. But, new libraries added to the 
> latest version such as kwant.kpm.conductivity or importing kwant.qsymm are 
> not working when I call them. What should I do?
>
> Best regards,
>
> Ali

Re: [Kwant] NNN hopping

[Anton Akhmerov]

Hi Yuhao,

The function "uniform" that you're using produces random numbers based
on its input. Therefore if you want these to be the same for two sites
in the same unit cell, you need to pass the unit cell coordinates as
input. These are available as site.tag . So combining everything
together you need

def nnn(site1, site2): return 1j *  (1-2*uniform(repr(site1.tag),salt))

Best,
Anton

On Fri, Mar 1, 2019 at 9:05 AM yuhao kang  wrote:
>
> Hello,
>
> I'm trying to add disorder in the NNN hopping term at the Haldane model.
>
> a, b = graphene.sublattices
> nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a))
> nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b))
> nnn_hoppings = nnn_hoppings_a + nnn_hoppings_b
>
> def nnn(site1, site2):
> return 1j *  (1-2*uniform(repr(site1),salt))
>
> sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings]] = nnn
>
> I want to vary the strength of nnn hopping in different unit cell, and keep 
> them to be the same in the same unit cell.
> But here, even the site a/b is in one unitcell, the nnn term for a and b are 
> different, is there any way to solve it?
>
> Thank you in advance!
>
> Best,
> yuhao

Re: [Kwant] Number of modes calculation

[Anton Akhmerov]

Hi Ran,

> 1- Does KWANT include the number of modes, M(E), inside the transmission 
> function or not. In other words does “smatrix.transmission(1,0)” contains 
> only T(E) or M(E)*T(E) ?

smatrix.transmission is indeed the sum of all the transmission
coefficients, so it is correctly proportional to the number of modes.
The conductance is equal to smatrix.transmission times the conductance
quantum.

Best,
Anton

Re: [Kwant] Problems of Running kwant.ldos(sys, 0.2) and kwant.plotter.map(sys, ldos)

[Anton Akhmerov]

Dear Zhanzhi,

I believe you missed the line "sys = sys.finalized()". Please use the
Kwant tutorial to learn Kwant instead, it's a much more systematic
source.

Best,
Anton

On Tue, Jan 1, 2019 at 11:13 PM Zhanzhi Jiang  wrote:
>
> Hi,
>
> I am following the vedio tutorial provided by Christoph  
> (https://kwant-project.org/doc/kwant-screencast-2014)
>
> When I run
>
> ldos = kwant.ldos(sys, 0.2)
>
> There is an error "TypeError: Expecting an instance of System"
>
> Then when I run
>
> kwant.plotter.map(sys, ldos)
>
> There is an error "NameError: name 'ldos' is not defined"
>
> So I cannot get the map of the local DOS as shown in the vedio. It looks like 
> the second error is a result of the first one. How could I solve this?
>
> Best,
> Zhanzhi
>

Re: [Kwant] Error visualizing the zigzag graphene lattice using kwant.plot() in tutorial

[Anton Akhmerov]

Hi Srilok,

Right now plotting of detached finalized leads isn't supported and is
tracked in the issue
https://gitlab.kwant-project.org/kwant/kwant/issues/18

Cheers,
Anton
On Tue, Nov 20, 2018 at 9:22 AM Srilok Srinivasan  wrote:
>
> Sorry I forgot to mention that I am using the Kwant version 1.4.01a1
>
> Thanks
>
> On Mon, Nov 19, 2018 at 5:38 PM Srilok Srinivasan  wrote:
>>
>>
>> Hello Kwant users,
>>
>> I am very new to kwant and I am still learning to use the code. So my 
>> apologizes for any trivial questions.
>>
>> I working on the tutorial notebooks I found in the here 
>> (https://kwant-project.org/mm16). In the tutorial 3.4.graphene_qshe.ipynb I 
>> am trying to visualize the zigzag graphene lattice created using the 
>> kwant.Builder and .shape() module using kwant.plot(sys). I am getting the 
>> following error
>>
>> zigzag_ribbon = kwant.Builder(kwant.TranslationalSymmetry([1, 0]))
>> #zigzag_ribbon[graphene.shape((lambda pos: abs(pos[1]) < 9 and pos[1]%2 
>> ==1), (0, 0))] = 0.2
>> #zigzag_ribbon[graphene.shape((lambda pos: abs(pos[1]) < 9 and pos[1]%2 
>> ==0), (0, 0))] = -0.2
>> zigzag_ribbon[graphene.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = 0
>> zigzag_ribbon[graphene.neighbors(1)] = 1
>> zigzag_ribbon.finalized()
>> kwant.plotter.plot(zigzag_ribbon)
>> #kwant.plotter.bands(zigzag_ribbon.finalized());
>>
>> ---
>> KeyError  Traceback (most recent call last)
>>  in 
>>   5 zigzag_ribbon[graphene.neighbors(1)] = 1
>>   6 zigzag_ribbon.finalized()
>> > 7 kwant.plotter.plot(zigzag_ribbon)
>>   8 #kwant.plotter.bands(zigzag_ribbon.finalized());
>>
>> ~/anaconda3/lib/python3.6/site-packages/kwant/plotter.py in plot(sys, 
>> num_lead_cells, unit, site_symbol, site_size, site_color, site_edgecolor, 
>> site_lw, hop_color, hop_lw, lead_site_symbol, lead_site_size, lead_color, 
>> lead_site_edgecolor, lead_site_lw, lead_hop_lw, pos_transform, cmap, 
>> colorbar, file, show, dpi, fig_size, ax)
>> 962
>> 963 site_size = make_proper_site_spec(site_size, fancy_indexing)
>> --> 964 site_color = make_proper_site_spec(site_color, fancy_indexing)
>> 965 site_edgecolor = make_proper_site_spec(site_edgecolor, 
>> fancy_indexing)
>> 966 site_lw = make_proper_site_spec(site_lw, fancy_indexing)
>>
>> ~/anaconda3/lib/python3.6/site-packages/kwant/plotter.py in 
>> make_proper_site_spec(spec, fancy_indexing)
>> 907 def make_proper_site_spec(spec, fancy_indexing=False):
>> 908 if callable(spec):
>> --> 909 spec = [spec(i[0]) for i in sites if i[1] is None]
>> 910 if (fancy_indexing and _p.isarray(spec)
>> 911 and not isinstance(spec, np.ndarray)):
>>
>> ~/anaconda3/lib/python3.6/site-packages/kwant/plotter.py in (.0)
>> 907 def make_proper_site_spec(spec, fancy_indexing=False):
>> 908 if callable(spec):
>> --> 909 spec = [spec(i[0]) for i in sites if i[1] is None]
>> 910 if (fancy_indexing and _p.isarray(spec)
>> 911 and not isinstance(spec, np.ndarray)):
>>
>> ~/anaconda3/lib/python3.6/site-packages/kwant/plotter.py in site_color(site)
>> 956 color_mapping = dict(zip(families, cycle))
>> 957 def site_color(site):
>> --> 958 return color_mapping[site.family]
>> 959 else:
>> 960 # Unknown finalized system, no sites access.
>>
>> KeyError: kwant.lattice.Monatomic([[1.0, 0.0], [0.5, 0.8660254037844386]], 
>> [0.0, 0.5773502691896258], '1', None)
>>
>> However, I am able to plot the dispersion and all that. Also, I am able to 
>> visualize the square lattice presented in other tutorials. Do we have to do 
>> something different to visualize non-square lattice ?
>>
>> Thanks for your help
>> --
>> Srilok Srinivasan
>> Graduate Student
>> Mechanical Engineering
>> Iowa State University, Ames, IA
>>
>
>
> --
> Srilok Srinivasan
> Graduate Student
> Mechanical Engineering
> Iowa State University, Ames, IA
>

Re: [Kwant] issues with Current calculation

[Anton Akhmerov]

On Fri, Oct 19, 2018 at 11:39 AM Christoph Groth  wrote:
>
> Sergey Slizovskiy wrote:

> > sudo add-apt-repository ppa:kwant-project/ppa
> > sudo apt-get update
> > sudo apt-get install kwant

I believe the package name is python3-kwant as written in
https://kwant-project.org/install#ubuntu-and-derivatives

> Kwant 1.3.3 and Tinyarray packages are available for Ubuntu bionic:
> https://launchpad.net/~kwant-project/+archive/ubuntu/ppa/+packages
>
> No idea what's going wrong with your installation of Linux Mint 19.
>
> Christoph

Re: [Kwant] Merging two systems

[Anton Akhmerov]

Dear Shivang Aragwal,

If the two systems do not contain same sites, then combining them is as
simple as syst1.update(syst2) (before finalizing).

See also
https://kwant-project.org/doc/1/reference/generated/kwant.builder.Builder#kwant.builder.Builder.update

Best,
Anton

On Wed, Oct 17, 2018 at 10:16 AM Shivang Agarwal <
shivang.agar...@iitgn.ac.in> wrote:

> Hello All,
>
> I was wondering if there is a way in Kwant to combine two systems. For
> example, lets say I build two systems sys1 and sys2 (on the kwant.plot
> graph, the two systems do not overlap, i.e. they are separated by a finite
> distance). Then, would it be possible to plot the two systems on the same
> graph, and then calculate transmission of combined sys1+sys2?
>
> Any help would be much appreciated.
>
> Best Regards,
> Shivang Agarwal
>
> --
> *Shivang Agarwal*
> Senior Undergraduate
> Coordinator - Academic Discussion Hours
> Discipline of Electrical Engineering
> IIT Gandhinagar
>
> Contact: +91-9869321451
>

Re: [Kwant] Bilayer material tight-binding in Kwant

[Anton Akhmerov]

Hi Marc,

Working in 2D is just fine in that case. I think it's also easier, so
unless you have a special need for 3D, I'd recommend to omit the
z-coordinate.

Best,
Anton
On Wed, Oct 10, 2018 at 12:44 PM Marc Vila  wrote:
>
> Dear Kwant developers,
>
>
> If I want to model a bilayer-type structure like bilayer graphene or 
> phosphorene, do I need to work in three dimensions being the perpendicular 
> lattice constant the interlayer spacing? Or can I still work in 
> two-dimensions since these materials are two-dimensional when the material is 
> periodic in the plane?
>
>
> Thanks in advance for the help!
>
>
> Kind regards,
>
>
> Marc

Re: [Kwant] seeking help

[Anton Akhmerov]

Dear Shyam Lochan Bora,

Right now your question is rather poorly formulated, and not suitable
for this mailing list. Please read this instruction on how to ask good
questions: https://stackoverflow.com/help/how-to-ask
Additionally, it appears that your problem is due the lack of the
Python knowledge, I recommend you to follow a Python course, or
otherwise learn the language more systematically.

Best regards,
Anton Akhmerov
On Fri, Aug 31, 2018 at 5:08 PM shyam lochan bora
 wrote:
>
> Dear sir,
> Please help me with the following error in the program
> from __future__ import print_function
> from ipywidgets import interact, interactive, fixed, interact_manual
> import ipywidgets as widgets
> import kwant
> from matplotlib import pyplot
> import numpy as np
> from holoviews.core.options import Cycle
> from types import SimpleNamespace
>
> def nanowire_chain():
> lat = kwant.lattice.chain()
> sys = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
>
> def onsite(onsite, p):
> return (2 * p.t - p.mu) * pauli.szs0 + p.B * 
> (np.cos(p.phi)*pauli.s0sz+ 
> np.sin(p.phi)*np.sin(p.theta)*pauli.s0sy+np.sin(p.phi)*np.sin(p.theta)*pauli.s0sx)+
>  p.delta * pauli.sxs0
>
> sys[lat(0)] = onsite
>
> def hop(site1, site2, p):
> return -p.t * pauli.szs0 - .5j * p.alpha * pauli.szsx
>
> sys[kwant.HoppingKind((1,), lat)] = hop
>
> return sys
>
>
> def spinful_kitaev_chain():
> lat = kwant.lattice.chain()
> sys = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
>
> def onsite(site, p):
> return (2 * p.t - p.mu) * pauli.szs0 + p.B * pauli.szsz
>
> sys[lat(0)] = onsite
>
> def hop(site1, site2, p):
> return -p.t * pauli.szs0 - 1j * p.delta * pauli.sys0
>
> sys[kwant.HoppingKind((1,), lat)] = hop
>
> return sys
>
>
> def find_gap(sys, p, resolution=1e-4):
> """Find gap in a system by doing a binary search in energy."""
>
> # This tells us if there are modes at a certain energy.
> if len(sys.modes(energy=0, args=[p])[0].momenta):
> return 0
>
> gap = step = min(abs(kwant.physics.Bands(sys, args=[p])(k=0))) / 2
> while step > resolution:
> step /= 2
> if len(sys.modes(gap, args=[p])[0].momenta):
> gap -= step
> else:
> gap += step
>
> return gap
>
>
> def spinorbit_band_gap(sys, mu, t, delta, Bs):
> sys = sys.finalized()
> alphas = [0.0, 0.1, 0.2, 0.3]
> p = SimpleNamespace(mu=mu, t=t, delta=delta,theta=theta,phi=phi)
>
> def gap(sys, p, alpha, B):
> p.alpha = alpha
> p.B = B
> return find_gap(sys, p)
>
> gaps = [gap(sys, p, alpha, B) for alpha in alphas for B in Bs]
> gaps = np.reshape(gaps, (len(alphas), -1))
> dims = {'kdims': [r'$B$'], 'vdims': ['Band gap']}
> B_crit = holoviews.VLine(np.sqrt(p.delta**2 + p.mu**2))
> plot = [holoviews.Curve((Bs, gaps[i]), label=r'$\alpha={}$'.format(
> alphas[i]), **dims) * B_crit for i, alpha in enumerate(alphas)]
> title = r'$\Delta={delta}$, $\mu={mu}$'.format(delta=np.round(p.delta, 
> 2), mu=np.round(p.mu, 2))
> style = {'xticks': [0, 0.1, 0.2, 0.3], 'yticks': [0, 0.05, 0.1], 
> 'fig_size': 150}
> plot = holoviews.Overlay(plot)
> return plot(plot=style)
>
>
> def title(p):
> try:
> title = r"$\alpha={alpha}$, $\mu={mu}$, $B={B}$, $\Delta={delta}$"
> title = title.format(alpha=np.round(p.alpha, 2),
>  mu=np.round(p.mu, 2),
>  B=np.round(p.B, 2),
>  delta=np.round(p.delta, 2))
> except AttributeError:
> title = r"$\mu={mu}$, $B={B}$, $\Delta={delta}$"
> title = title.format(mu=np.round(p.mu, 2),
>  B=np.round(p.B, 2),
>  delta=np.round(p.delta, 2))
> return title
>
> style = {'k_x': np.linspace(-1, 1, 101),
>  'xdim': r'$k$',
>  'ydim': r'$E/t$',
>  'xticks': [-1, 0, 1],
>  'yticks': [-1, 0, 1],
>  'xlims': [-1, 1],
>  'ylims': [-1.5, 1.5],
>  'title': title}
> sys = nanowire_chain()
> p = SimpleNamespace(t=1, mu=0.1, delta=0.1, B=0.3,theta=np.pi/2,phi=np.pi/2)
> alphas = np.linspace(0, 0.4, 10)
> holoviews.HoloMap({alpha: spectrum(sys, p.update(alpha=alpha), **style) for 
> alpha in alphas}, kdims=[r'$\alpha$'])
>
> c:\python36\lib\site-packages\kwant\linalg\lll.py:103: FutureWarning: `rcond` 
> parameter will change to the default of machine precision times ``max(M, N)`` 
&

Re: [Kwant] next nearest neighbour

[Anton Akhmerov]

Dear KS Chan,

You can use kwant.wraparound.wraparound [1] and kwant.plotter.spectrum
[2] to compute and plot the band structure of a lead with longer than
nearest neighbor hoppings.

When attaching such a lead to the scattering region in the latest
version Kwant will automatically double the unit cell, which
unfortunately means that you need to create two different leads: one
for band structure and another one for scattering.

Best,
Anton

[1]: 
https://kwant-project.org/doc/1/reference/generated/kwant.wraparound.wraparound#kwant.wraparound.wraparound
[2]: 
https://kwant-project.org/doc/1/reference/generated/kwant.plotter.spectrum#kwant.plotter.spectrum
On Tue, Aug 14, 2018 at 6:19 PM Prof. CHAN Kwok Sum
 wrote:
>
> Dear Kwant Developer,
>
> I note there was a discussion in the community about next nearest neighbor 
> hopping in the community (Leads with greater than NN hopping).
>
> However, I still have the following questions about next nearest neighbour 
> hopping in the lead after reading the discussion.
>
> 1.Can we add next nearest neighbour hopping to a lead without expanding 
> the unit cell, if there are only hoppings between nearest neighbouring unit 
> cells even with next nearest neighbour hopping? Can the band structure be 
> calculated properly using the Kwant library?
>
> 2.When the lead is attached to the system, will the next nearest 
> neighbour hopping be included in the action automatically?
>
> Thank you in advance for your help.
>
> 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.

Re: [Kwant] Calculating valley polarized conductances using Kwant

[Anton Akhmerov]

Hi Kevin,

> I just wanted to verify that the indices of the modes I am finding match up 
> with the indices of the submatrix.

This is indeed correct. The sort order of the scattering matrix
elements matches the sort order of the modes in the lead_info.

Best,
Anton

Re: [Kwant]

[Anton Akhmerov]

Dear KS Chan,

> When we use plot to plot a 3D system, is it possible to remove, the x, y, z 
> axes in the final figure? In 2D, no x, y axes is shown. It is nice to have 
> this flexibility.

This can be done by using regular matplotlib tools, for example like this:

fig = kwant.plot(syst, ..., show=False)
fig.axes[0].set_axis_off()
pyplot.show()

Re: [Kwant] Finding hamiltonian of system (error in user defined function 'onsite')

[Anton Akhmerov]

Dear Shivang,

You must supply to the hamiltonian_submatrix extra arguments that
onsite expects. Inspecting the error message in more detail would
reveal that. In order to figure out how to provide the extra arguments
to hamiltonian_submatrix, please consult the documentation.

Best,
Anton
On Thu, Jul 5, 2018 at 7:18 PM Shivang Agarwal
 wrote:
>
> Hello All,
>
> I want to obtain the hamiltonian for my system, for which I am using the 
> sys.hamiltonian_submatrix() command, but I get the following error.
>
> UserCodeError: Error occurred in user-supplied value function "onsite".
>
> Below is the code that I am using.
>
> lat = kwant.lattice.honeycomb(a = a)
> p, q = lat.sublattices
>
> def make_system():
>
> def hopping(sitei, sitej, phi, salt):
> xi, yi = sitei.pos
> xj, yj = sitej.pos
> return -t * np.exp(-0.5j * phi * (xi - xj) * (yi + yj))
>
> def onsite(site, phi, salt):
> x, y = site.pos
> return 0
>
> def central_region(pos):
> x, y = pos
> return abs(x) < 2*a and abs(y) < 3*a
>
> sys = kwant.Builder()
> sys[lat.shape(central_region, (0, 0))] = onsite
> sys[lat.neighbors()] = hopping
>
> sym = kwant.TranslationalSymmetry((-a, 0))
> lead = kwant.Builder(sym)
> lead[lat.shape(lambda s: abs(s[1]) < 2*a, (0, 0))] = 0
> lead[lat.neighbors()] = hopping
>
> sys.attach_lead(lead)
> sys.attach_lead(lead.reversed())
>
> return sys.finalized()
>
> sys = make_system()
> sparse_mat = sys.hamiltonian_submatrix()  > Throws the error
>
> What confuses me is that the code works fine (when I do further calculations 
> on sys) if I comment the hamitonian_submatrix() command but throws up the 
> error otherwise. That made me think that the code itself should be correct, 
> but I could not figure out why that line gives an error.
>
> I would be glad if someone can help me out with the same. Any help will be 
> much appreciated!
>
> Sincerely,
> Shivang Agarwal
>
> --
> Shivang Agarwal
> Junior Undergraduate
> Discipline of Electrical Engineering
> IIT Gandhinagar
>
> Contact: +91-9869321451

Re: [Kwant] AttributeError: 'FiniteSystem' object has no attribute 'site' [site = sys.site(i) not working]

[Anton Akhmerov]

Dear Shivang,

I believe the code you consider is from a previous version of Kwant.
Now it'd be syst.sites[i]. Can you please let us know where you found
that code? Was it in the Kwant documentation?

Best,
Anton
On Tue, Jun 12, 2018 at 8:58 PM Shivang Agarwal
 wrote:
>
> Hello All,
>
> I am looking to plot the wavefunction (both magnitude and phase on separate 
> plots).
> The kwant.plotter.map method doesn't give a very appealing look as it plots 
> only in large sized triangles. The kwant.plotter.plot method seems better but 
> I get stuck in one part. There is a piece of code that reads:
>
> def family_shape(i):
> site = sys.site(i)
> return ('p', 3, 180) if site.family == a else ('p', 3, 0)
>
> However, it gives the following error:
>
> AttributeError: 'FiniteSystem' object has no attribute 'site'
>
> This snippet was taken from online kwant tutorials itself. I am unable to get 
> to the bottom of this error message because the tutorial uses the exact same 
> code. Is there an update to some library of kwant that solves this? Or is 
> there any other way to get better plots for wavefunctions on a system without 
> using the sys.site(i) command?
> I would like the wave functions to look smooth and continuous.
>
> Any help would be appreciated.
>
> Best Regards,
> Shivang
>
> --
> Shivang Agarwal
> Junior Undergraduate
> Discipline of Electrical Engineering
> IIT Gandhinagar
>
> Contact: +91-9869321451

Re: [Kwant] Reflectionless contacts

[Anton Akhmerov]

Hi Eleni,

It's a hard problem, and usually we reduce reflection by creating a
smooth interface between the scattering region and the contact. As far
as measuring the quality of such interface goes, I think the most
practical approach is to replace the scattering region with another
lead with dispersion mimicking that of the scattering region as
closely as possible.

Best,
Anton

On Sun, Apr 22, 2018 at 12:13 PM,   wrote:
> Hello everyone,
>
> I have defined two semi-infinite contacts in my system as follows:
>
> 
> def lead_shape(p):
>   x, y, z = p
>   return -10 < x < 10 and -6 < y < 6 and 0 < z < 30
>
> sym_lead = kwant.TranslationalSymmetry(lattice.vec((-6, 0, 0)))
> lead = kwant.Builder(sym_lead)
> lead[lattice.shape(lead_shape, (0, 0, 15))] = 0
> 
>
> Is there an easy way to calculate any reflections from the lead back into
> the conductor?
> I mean, ideally, I would like all electrons entering the contacts to stay
> there.
>
>
> Regards
>
>
> Eleni
>
>
> --
> Dr. Eleni Chatzikyriakou
> Computational Physics lab
> Aristotle University of Thessaloniki
> elch...@auth.gr - tel:+30 2310 998109
>

Re: [Kwant] Visualize current through a 2D cut for a 3D system

[Anton Akhmerov]

Dear Eleni,

I think this should be doable with not too much of low level work.
* Firstly you probably will want to select only the hoppings that are
close to the 2D cut. You can do this by using "where" parameter to
Current operator (see
https://kwant-project.org/doc/1/reference/generated/kwant.operator.Current#kwant.operator.Current)
* After computing the current across those hoppings, you would need to
call kwant.plotter.interpolate_current. This function only does the
interpolation, but no plotting. Having that interpolation in 3D you
need to slice it and probably take a projection onto 2D plane.
* Finally, the resulting 2D array you can feed to kwant.plotter.streamplot

kwant.plotter.current is essentially a straightforward combination of
kwant.plotter.interpolate_current and kwant.plotter.streamplot

Also: your problem sounds like a useful thing to do. Please let me
know if you succeed, and what you ended up doing (or if you face any
further problems as well).

Best,
Anton

On Fri, Mar 23, 2018 at 5:22 PM,   wrote:
> Hello everyone,
>
> Is there currently way to visualize the current through a 2D cut for a 3D
> system that does not involve messing with the kwant code?
>
> I have a 3D system and would like to see the differences in the current when
> I switch some specific hoppings on and off.
>
>
> Regards,
>
>
>
> --
> Dr. Eleni Chatzikyriakou
> Computational Physics lab
> Aristotle University of Thessaloniki
> elch...@auth.gr - tel:+30 2310 998109
>

Re: [Kwant] combine two kwant.plotter.current figures into a single one

[Anton Akhmerov]

Dear Tibor,

All Kwant plotting functions accept an Axes instance as an input
argument, in which case they will modify that axes. Therefore to
achieve what you want you need to prepare a figure with multiple axes
and pass these axes to the Kwant plotting functions as input.

Best,
Anton

On Tue, Jan 30, 2018 at 4:26 PM, Tibor Sekera  wrote:
> Dear all,
>
> sorry, this might be more a matplotlib question than a kwant question.
>
> Below is a minimal code which produces two figures, each showing a different
> dataset: a current in the scattering region due to a wavefunction in a lead.
>
> I would like to combine the two into a single plot, distinguishing the two
> datasets by e.g. different colormaps (and two different colorbars).
>
> Any ideas?
>
> Thanks a lot,
> Tibor
>
> __
>
> import kwant
> from math import pi, sqrt
> from cmath import exp
> import numpy
> from matplotlib import pyplot as plt
>
> #some parameters
> t=1
> L=30
> W=30
> EF=0.1
> phi=0.01
>
> def hopping(sitei, sitej, phi):
> xi, yi = sitei.pos
> xj, yj = sitej.pos
> return -t*exp(1j*2*pi*phi*2/sqrt(3)*(yj - yi)*(xi + xj)/2)
>
> def make_system():
>
> lat = kwant.lattice.general([(sqrt(3)*1/2, 1/2), (0, 1)],
>  [(0, 0), (1/(2*sqrt(3)),1/2)], norbs=1)
>
> def central_region(pos):
> x, y = pos
> return abs(x) < L/2 and abs(y) < W/2
>
> def lead0_shape(pos):
> x, y = pos
> return abs(x) < L/2
>
> #scattering region
> sys = kwant.Builder()
> sys[lat.shape(central_region, (0,0))] = -EF
> sys[lat.neighbors()] = hopping
> sys.eradicate_dangling()
>
> #leads
> sym = kwant.TranslationalSymmetry(lat.vec((0,1)))
> lead0 = kwant.Builder(sym)
> lead0[lat.shape(lead0_shape, (0,0))] = -EF
> lead0[lat.neighbors()] = hopping
> lead0.eradicate_dangling()
> lead1=lead0.reversed()
>
> return sys, lead0, lead1
>
> sys, lead0, lead1 = make_system()
> sys.attach_lead(lead0)
> sys.attach_lead(lead1)
> sysf = sys.finalized()
>
>
> wf = kwant.wave_function(sysf, energy=0.0001, params=dict(phi=phi))
> J = kwant.operator.Current(sysf)
>
> psi1 = wf(0)[0]
> e_current1 = J(psi1, params=dict(phi=phi))
>
> psi2 = wf(1)[0]
> e_current2 = J(psi2, params=dict(phi=phi))
>
> kwant.plotter.current(sysf, e_current1, cmap="OrRd", colorbar=True,
> show=False)
> ax = plt.gca()
> plt.colorbar(ax.get_children()[-2], ax=ax)
>
> kwant.plotter.current(sysf, e_current2, cmap="Blues", colorbar=True,
> show=False)
> ax = plt.gca()
> plt.colorbar(ax.get_children()[-2], ax=ax)
>
> plt.show()

Re: [Kwant] getting less number of peaks in transmission for 1D chain

[Anton Akhmerov]

Hi Sudin,

Could it be that some of the peaks are overlapping because there are
degeneracies in the spectrum?

Best,
Anton

On Tue, Jan 23, 2018 at 4:53 AM, SUDIN GANGULY  wrote:
> Dear Sir,
>
> Recently I was doing a problem on a 1D chain. The problem I faced is as
> follows.
>
> For a 1D chain, suppose I have taken 5 sites and set the hopping strength
> between lead and system very small compared to the system hopping.  Now
> since we have 5 sites, from the Hamiltonian of the system we'll have 5
> eigenvalues and if we calculate the transmission coefficient we shall have 5
> peaks each occurs at particular eigenvalue.
>
> But in my case, I am getting only three peaks. Am I missing something here?
> Please help me. The following code I have used for the transmission
> calculation.
>
> #===CODE==
> import kwant
> # For plotting
> from matplotlib import pyplot
> import numpy as np
>
>
> def make_system(a=1, t=1.0, L=5):
> # Start with an empty tight-binding system and a single square lattice.
> # `a` is the lattice constant (by default set to 1 for simplicity.
>
> lat = kwant.lattice.chain(a)
>
> syst = kwant.Builder()
>
>  Define the scattering region. 
> for x in range (L):
> syst[(lat(x))] = 0
>
> #syst[(lat(x) for x in range(L))] = onsite
> syst[lat.neighbors()] = -t
>
>
>  Define and attach the leads. 
> lead = kwant.Builder(kwant.TranslationalSymmetry([-a]))
> lead[(lat(0))] = 0
> lead[lat.neighbors()] = -t
>
> syst.attach_lead(lead)
> syst.attach_lead(lead.reversed())
>
> #== set hopping between lead and system 
> for i in syst.leads[0].interface:
> for j in syst.neighbors(i):
> syst[i,j] = 0.1
>
> for i in syst.leads[1].interface:
> for j in syst.neighbors(i):
> syst[i,j] = 0.1
>
> return syst
> #
> def main():
> syst = make_system()
>
> # Check that the system looks as intended.
> kwant.plot(syst)
>
> # Finalize the system.
> syst = syst.finalized()
>
>
> Es = np.linspace(-1.99, 1.99, 800)
>
> data = []
> totdos = []
> for e in Es:
> smatrix = kwant.smatrix(syst, e)
> T = smatrix.transmission(1, 0)
> ldos = kwant.ldos(syst, e)
> totdos.append(sum(ldos)/len(ldos))
> data.append(T)
>
> print ('%0.3f' % e, '%.4f' % T, sum(ldos)/len(ldos))
>
> pyplot.figure()
> pyplot.plot(Es, data)
> pyplot.xlabel("energy [t]")
> pyplot.ylabel("conductance [e^2/h]")
> pyplot.show()
>
>
> # Call the main function if the script gets executed (as opposed to
> imported).
> # See .
> if __name__ == '__main__':
> main()
>
>
> #End CODE===
> --
> সুদিন

Re: [Kwant] honeycomb lattice scattering region with square lattice leads

[Anton Akhmerov]

Dear Hosein Khani,

In Kwant the hoppings between the scattering region and the leads are
the same as the lead hoppings. I hope this answers your question.

Best,
Anton

On Mon, Dec 11, 2017 at 12:45 AM, Khani Hosein  wrote:
> Dear developers,
> I am considering a system with honeycomb lattices (scattering region)
> connected to square lattices (leads). For the square lattices of the leads,
> I want to use the honeycomb lattice but with hoppings same as square
> lattices, so it works like square lattice leads. My question is: how can I
> find out the hoppings between the scattering region and the leads. The leads
> are attached to the scattering region using sys.attach_lead(lead).
> Regards,
> Hosein Khani
>

Re: [Kwant] Zero energy states with KPM

[Anton Akhmerov]

Dear Antonio,

KPM by itself doesn't allow to adjust the resolution locally. If you
need detailed information about a small part of the spectrum, sparse
diagonalization is a more appropriate tool.

Best,
Anton

On Tue, Nov 28, 2017 at 1:58 PM, Antonio Lucas Rigotti Manesco
 wrote:
> Dear Kwant users,
>
> I am trying to obtain the local density of states of a topological (finite)
> system by using KPM approximation. My problem is that the gap is too small
> (as I am considering realistic parameters) and this results in more states
> than the topologically protected ones appearing for zero energy LDoS. If I
> control the energy resolution the memory consumption increases too much.
>
> So my question is if there is any parameter I can control to have a better
> energy resolution only near zero energy or maybe somewhere I can compensate
> the memory consumption (for example, losing resolution at the local
> distribution in order to obtain better energy resolution).
>
> Thanks in advance.
> --
> Antônio Lucas Rigotti Manesco
> PhD fellow - University of São Paulo, Brazil

Re: [Kwant] Direction of wavevector of kwant.physics.Bands?

[Anton Akhmerov]

Hi George,

We (the Kwant developers) realized that the mode ordering isn't very
clear, and are tracking the fixes in this issue:
https://gitlab.kwant-project.org/kwant/kwant/issues/143

Please take a look at the related changes in the documentation. If you
find those not sufficiently clear, please let us know (e.g. by joining
the issue discussion or replying here).

Best,
Anton

On Fri, Nov 10, 2017 at 4:37 PM, George Datseris
 wrote:
> Hello,
>
> Let's say I have a system, and a lead attached to the left side of the
> system extending to infinity towards the left.
>
> As far as kwant is concerned, the translation invariance vector for the lead
> points towards the left, i.e. in negative x direction.
>
> Now I call `kwant.physics.Bands` or equivalently `kwant.plotter.bands` on
> this lead.
>
> I get a plot of the energy bands, with the wavevector ranging from -π to +π
> , in units of inverse lattice constants of the lead.
>
> My question is the following: What does the positive x-axis of this plot
> correspond to? Does positive wavevector in this plot correspond to positive
> x-direction (opposite to the lead translational vector) or to negative
> x-direction, in parallel with the lead translational vector?
>
> I know that e.g. incoming modes must have negative group velocity, but this
> doesn't answer my question, since I am using a graphene zigzag lead, which
> means that I have negative group velocities for both positive and negative
> wavevectors.
>
> Also, I have checked the following documentations trying to answer my
> question:
>
> kwant.physics.Bands
> PropagatingModes
> FAQ of development documentation page
> kwant.plotter.band
> kwant.physics.modes
>
> It is still unclear to me whether in my case, positive wavevector
> corresponds to the positive x-direction of my system or the negative
> x-direction of my system, due to the lead convention.
>
> If this is stated clearly somewhere in the kwant documentation, please
> excuse me for not finding it, and be kind enough to point me there.
>
> P.S.: I need this because I am trying to calculate from which graphene
> "valley" each incoming mode belongs to. If the k_x wavevector of the
> incoming valley is postive, then it has valley pseudo-spin ξ=1, but I have
> to know whether my "positive wavevector" and the lead's positive wavevectors
> are the same or opposite.
>
> Best regards,
> George Datseris

Re: [Kwant] Majorana polarization

[Anton Akhmerov]

Hi Antonio.

I imagine you could evaluate that vector field at the positions of the
lattice sites and then feed it to an interpolation routine.

Best,
Anton

On Tue, Oct 31, 2017 at 12:35 PM, Antonio Lucas Rigotti Manesco
<antonio...@usp.br> wrote:
> Hi,
>
> As Majorana polarization is a vector field, I am trying to plot the values
> using quiver by adding the two functions:
>
> def plot_vector_field(sys, Maj_x, Maj_y):
> xmin, ymin = min(s.tag for s in sys.sites)
> xmax, ymax = max(s.tag for s in sys.sites)
> x, y = np.meshgrid(np.arange(xmin, xmax+1), np.arange(ymin, ymax+1))
>
> plt.quiver(x, y, Maj_x, Maj_y)
> plt.show()
>
> def majorana_polarization(sys, energy):
> wf = kwant.wave_function(sys, energy)
> psi = wf(0)
> eup, edown, hup, hdown = psi[::2], psi[1::2], psi[3::2], - psi[2::2]
> Maj = eup * hup + edown * hdown
> print(Maj)
> Maj_x = np.real(Maj)
> Maj_y = np.imag(Maj)
> plot_vector_field(sys, Maj_x, Maj_y)
>
> But since I am dealing with a honeycomb lattice the meshgrid will not work.
> How could I generate that using the sys sites?
>
> 2017-10-30 19:10 GMT-02:00 Anton Akhmerov <anton.akhmerov...@gmail.com>:
>>
>> Dear Antonio,
>>
>> In Kwant 1.3 you can compute expectation values of any local operator,
>> applied to any wave function, as described in this Kwant tutorial:
>> https://kwant-project.org/doc/1/tutorial/operators
>>
>> It is also possible that you can compute what you want using the
>> output of ldos. The ordering of internal degrees of freedom of each
>> there is the same as in the Hamiltonian you define.
>>
>> Best,
>> Anton
>>
>> On Mon, Oct 30, 2017 at 9:44 PM, Antonio Lucas Rigotti Manesco
>> <antonio...@usp.br> wrote:
>> > Dear Kwant users,
>> >
>> > I am interested to calculate Majorana polarization for a system I am
>> > currently studying. I was wondering if it is possible to obtain it from
>> > kwant.ldos. The system is constructed as a kronecker product between
>> > particle-hole  and spin in that order but I am not sure how kwant.ldos
>> > orders the internal degrees of freedom. Does any of you have some hint
>> > about
>> > that?
>> >
>> > Thanks in advance.
>
>
>
>
> --
> Antônio Lucas Rigotti Manesco
> PhD fellow - University of São Paulo, Brazil

Re: [Kwant] Majorana polarization

[Anton Akhmerov]

Dear Antonio,

In Kwant 1.3 you can compute expectation values of any local operator,
applied to any wave function, as described in this Kwant tutorial:
https://kwant-project.org/doc/1/tutorial/operators

It is also possible that you can compute what you want using the
output of ldos. The ordering of internal degrees of freedom of each
there is the same as in the Hamiltonian you define.

Best,
Anton

On Mon, Oct 30, 2017 at 9:44 PM, Antonio Lucas Rigotti Manesco
 wrote:
> Dear Kwant users,
>
> I am interested to calculate Majorana polarization for a system I am
> currently studying. I was wondering if it is possible to obtain it from
> kwant.ldos. The system is constructed as a kronecker product between
> particle-hole  and spin in that order but I am not sure how kwant.ldos
> orders the internal degrees of freedom. Does any of you have some hint about
> that?
>
> Thanks in advance.

Re: [Kwant] Connecting leads of a different lattice type.

[Anton Akhmerov]

On Fri, Oct 20, 2017 at 4:45 PM, Dripto Debroy <driptodeb...@gmail.com> wrote:
> That makes a lot of sense, would I be correct in assuming that I would need
> to define new lattices for each of the rows of square matrices before adding
> the hoppings by hand?

You only need to define two lattices: one honeycomb, and one square
for the lead.

Best,
Anton

>
> Much appreciated.
> Dripto
>
> On Oct 20, 2017 5:52 AM, "Anton Akhmerov" <anton.akhmerov...@gmail.com>
> wrote:
>>
>> Hi Dripto,
>>
>> The correct approach to implement that is to do the following:
>> 1. add a sufficient number of sites from the square lattice to the
>> scattering region, so that the lead may be attached to those (this
>> should amount to 1 unit cell of the lead).
>> 2. manually connect those sites to the honeycomb lattice. Since there
>> are too many possibilities to do so, there's no way to make this
>> choice automatically.
>> 3. attach the lead as usual using Builder.attach_lead.
>>
>> Best,
>> Anton
>>
>> On Thu, Oct 19, 2017 at 8:31 PM, Dripto Debroy <driptodeb...@gmail.com>
>> wrote:
>> > Hello all,
>> >
>> > I am hoping to use KWANT to model a graphene flake with two metal leads.
>> > The
>> > issue I am running in to is figuring out how to connect leads which have
>> > a
>> > square lattice and graphene which has a honeycomb lattice. Is this
>> > something
>> > that KWANT can handle? I was unable to find any documentation on mixing
>> > lattice types and modeling heterojunctions in general.
>> >
>> > Thanks for your help,
>> > Dripto

Re: [Kwant] Connecting leads of a different lattice type.

[Anton Akhmerov]

Hi Dripto,

The correct approach to implement that is to do the following:
1. add a sufficient number of sites from the square lattice to the
scattering region, so that the lead may be attached to those (this
should amount to 1 unit cell of the lead).
2. manually connect those sites to the honeycomb lattice. Since there
are too many possibilities to do so, there's no way to make this
choice automatically.
3. attach the lead as usual using Builder.attach_lead.

Best,
Anton

On Thu, Oct 19, 2017 at 8:31 PM, Dripto Debroy  wrote:
> Hello all,
>
> I am hoping to use KWANT to model a graphene flake with two metal leads. The
> issue I am running in to is figuring out how to connect leads which have a
> square lattice and graphene which has a honeycomb lattice. Is this something
> that KWANT can handle? I was unable to find any documentation on mixing
> lattice types and modeling heterojunctions in general.
>
> Thanks for your help,
> Dripto

Re: [Kwant] Importing Hamiltonian from file

[Anton Akhmerov]

Hi Marta,

I would like to know a bit more details about the problem you
encounter. Do you need to read a file with a specific format from
Python? What do you mean by the "full matrix", is it the sparse matrix
representation of the complete Hamiltonian, is it a list of onsite and
hopping values, or is it something entirely different?

Best,
Anton

On Mon, Oct 16, 2017 at 3:26 PM, Marta Brzezińska
 wrote:
> Good afternoon,
> starting from a tight-binding Hamiltonian saved in external file, how can I 
> import the full matrix to Kwant in order to perform transport calculations?
>
> Best regards,
> Marta Brzezinska

Re: [Kwant] 2D Current Density Streamplot of 3D System

[Anton Akhmerov]

Hi Mitchell,

Looks good, thanks for sharing. Are you satisfied with the usefulness of
such a plot?

By the way you could get a better looking result by using any plotting
library that supports true 3D graphics (matplotlib does not do actual 3D
rendering). Within Python world, possible options are plotly, mayavi, or
ipyvolume.

Cheers,
Anton

On Wed, Oct 4, 2017 at 10:23 AM, Mitchell Greenberg <
mitchell.greenb...@my.jcu.edu.au> wrote:

> Hi Anton,
>
>
> Had a go with getting the quiver plot working and from what I get it seems
> I've produced a graph that makes sense. My only thing I'm not 100% sure
> I've done right is the handling of the negative values; I assume that if I
> have site a and site b with a negative current value it means current is
> flowing from b to a, whereas if the current value is positive current is
> flowing from a to b.
>
>
> In order to account for this in the graph I plot each quiver point
> separately, see if the current is negative and if it is change the arrow to
> point into the site. I account for the magnitude of the current by setting
> the line width to each point to be proportional to the absolute value of
> the current. The attached code is below:
>
>
> #Plot current over hoppings
> wf = kwant.wave_function(syst, E,args=[gate])
> print('Calculated Wavefunctions')
> J_0 = kwant.operator.Current(syst)
> psi = wf(0)[0]
>
> current = J_0(psi)
>
> orderedcurrent = [(syst.pos(i), syst.pos(j), k) for (i, j), k in
> zip(syst.graph, current)]
>
> #From ordered current we know that the first array contains the
> start location of the hopping
> #and the second array contains the end location thus we can get
> XYZ,UVW from this.
>
> X = [orderedcurrent[i][0][0] for i in range(len(orderedcurrent))]
> Y = [orderedcurrent[i][0][1] for i in range(len(orderedcurrent))]
> Z = [orderedcurrent[i][0][2] for i in range(len(orderedcurrent))]
> U = [orderedcurrent[i][1][0] for i in range(len(orderedcurrent))]
> V = [orderedcurrent[i][1][1] for i in range(len(orderedcurrent))]
> W = [orderedcurrent[i][1][2] for i in range(len(orderedcurrent))]
>
>
> #Get the current magnitude at each hopping
> Mag = [orderedcurrent[i][2] for i in range(len(orderedcurrent))]
>
>
> fig = pyplot.figure(figsize=(12,10),dpi=100)
> ax = fig.gca(projection='3d')
> #Plot each point individually, if the current magnitude is
> negative, plot  the arrow head going into the site rather than out
> for i in range(len(orderedcurrent)):
> if Mag[i] < 0:
> ax.quiver(X[i], Y[i], Z[i], U[i], V[i], W[i], length =
> (abs(Mag[i])), pivot="tip")
> else:
> ax.quiver(X[i], Y[i], Z[i], U[i], V[i], W[i], length =
> (abs(Mag[i])), pivot="tail")
> pyplot.show()
>
>
> Thank You,
>
> Mitchell Greenberg
>
> --
> *From:* anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of
> Anton Akhmerov <anton.akhmerov...@gmail.com>
> *Sent:* Tuesday, October 3, 2017 11:22:26 PM
> *To:* Mitchell Greenberg
>
> *Cc:* kwant-discuss@kwant-project.org
> *Subject:* Re: [Kwant] 2D Current Density Streamplot of 3D System
>
> Hi Mitchell,
>
> There's nothing readily available, but it should certainly be doable.
> Firstly you need to calculate the coordinates of all sites, combine this
> with the graph information, and feed all this data to e.g. the matplotlib's
> quiver plot ( https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#
> mpl_toolkits.mplot3d.Axes3D.quiver ). See the recent reply of Joe (
> https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01425.html
> ) for assemblying the coordinates of all hoppings.
>
> Please share your code and your experience with this kind of plot if you
> succeed—do you think it would be universally useful?
>
> Best,
> Anton
>
> On Tue, Oct 3, 2017 at 3:15 PM, mitchell.greenberg <
> mitchell.greenb...@my.jcu.edu.au> wrote:
>
>> Hi Joseph and Anton,
>>
>> Today I had a talk with my supervisor and he wondered if it was possible
>> to plot the raw current data before interpolation, just to see how current
>> travels through each of the hoppings.
>>
>> From previous questions I know that the current array returned from J0
>> (psi) is in the same from as the system graph, however I was wondering if
>> there was a simple way to plot this array to see the current vectors at
>> each site position.
>>
>> Sorry if this question is naive.
>>
>> Thank 

Re: [Kwant] About the KPM method

[Anton Akhmerov]

Dear Kuangyia,

KPM method can be efficiently applied to any tight-binding system,
including infinite ones. Unfortunately, I didn't encounter an
implementation of KPM suitable for these geometries (e.g. it has to
deal with variable-sized vectors). In addition to that, the Kwant's
implementation of KPM is relatively basic, it does not implement
possible optimizations and it does not implement several useful
generalizations of KPM.

Best,
Anton

On Wed, Oct 4, 2017 at 10:53 AM, kuangyia lee  wrote:
> Dear Anton and other Kwant users,
>
>
> I am currently using the KPM method to calculate the density of states:
>
> spectrum = kwant.kpm.SpectralDensity(fsyst)
> dos = spectrum(energies)
>
>
> In the tutorial, this method is successfully applied to calculate the dos of
> a large closed system. But I am wondering that here it can also be applied
> to the case of an open system (i.e., a system composed of scattering region
> and leads).
>
> Thanks a lot !
>
> Best regards,
> Kuangyia

Re: [Kwant] 2D Current Density Streamplot of 3D System

[Anton Akhmerov]

Hi Mitchell,

There's nothing readily available, but it should certainly be doable.
Firstly you need to calculate the coordinates of all sites, combine this
with the graph information, and feed all this data to e.g. the matplotlib's
quiver plot (
https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#mpl_toolkits.mplot3d.Axes3D.quiver
). See the recent reply of Joe (
https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01425.html
) for assemblying the coordinates of all hoppings.

Please share your code and your experience with this kind of plot if you
succeed—do you think it would be universally useful?

Best,
Anton

On Tue, Oct 3, 2017 at 3:15 PM, mitchell.greenberg <
mitchell.greenb...@my.jcu.edu.au> wrote:

> Hi Joseph and Anton,
>
> Today I had a talk with my supervisor and he wondered if it was possible
> to plot the raw current data before interpolation, just to see how current
> travels through each of the hoppings.
>
> From previous questions I know that the current array returned from J0
> (psi) is in the same from as the system graph, however I was wondering if
> there was a simple way to plot this array to see the current vectors at
> each site position.
>
> Sorry if this question is naive.
>
> Thank You,
> Mitchell Greenberg
>
>
> Sent from Samsung Mobile
>
>
>  Original message 
> From: "mitchell.greenberg"
> Date:17/09/2017 7:32 PM (GMT+10:00)
> To: Anton Akhmerov
> Cc: kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>
> This sender failed our fraud detection checks and may not
> be who they appear to be. Learn about spoofing
> <http://aka.ms/LearnAboutSpoofing>
> Feedback <http://aka.ms/SafetyTipsFeedback>
> Hi Anton,
>
> Thanks for that. While my lattice is polyatomic I just switched site.pos
> to site.tag in the code I posted last message to get the tags for my
> lattice. Doing that I see the tags for the z axis range from 0 to 9 (makes
> sense as I have 10 layers), which explains why 5.3 worked when I used it
> but not 10. My function runs now an accurately deletes the sites without
> any problems.
>
> Thank You,
> Mitchell Greenberg
>
>
> Sent from Samsung Mobile
>
>
>  Original message 
> From: Anton Akhmerov
> Date:17/09/2017 5:11 AM (GMT+10:00)
> To: Mitchell Greenberg
> Cc: Joseph Weston ,kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>
> Hi Mitchel,
>
> Deleting and adding sites is done by their lattice coordinates
> (site.tag), not their real space position (site.pos); if you want to
> access a site close to a certain position, you can use the method
> Monatomic.closest (see
> https://kwant-project.org/doc/1/reference/generated/kwant.
> lattice.Monatomic#kwant.lattice.Monatomic.closest).
>
> Best,
> Anton
>
> On Sat, Sep 16, 2017 at 2:34 PM, Mitchell Greenberg
> <mitchell.greenb...@my.jcu.edu.au> wrote:
> > Hi Joseph and Anton,
> >
> >
> > Thank you for the advice with getting the 2d graph up and running, its
> now
> > producing results that look good. Sorry that I haven't replied in over 2
> > weeks, uni has gotten into full swing.
> >
> >
> > I wasn't sure if I needed to start a new topic for this (if I need to
> let me
> > know) but I am now trying to remove lattice points randomly from my
> sample.
> > I have created a function that can remove random lattice points (it is
> > rather messy) however I have noticed that when I try to delete sites for
> > layers other than at z=0 and z=5.3 I get the following error.:
> >
> >
> > KeyError: Site(kwant.lattice.Monatomic([[4.3763, 0.0, 0.0], [0.0,
> 3.3136,
> > 0.0], [0.0, 0.0, 5.3]], [0.0, 0.0, 0.0], '0', 1), array([50, 50, 10]))
> >
> >
> > This error seems to pop up when I try to delete a site that doesn't
> exist at
> > the specified coordinates however when I use:
> >
> > Positions=[syst.sites[i].pos for i in range(len(syst.sites))]
> >
> >
> > I see that sites definitely exist for z=10 (should be 10.6 but the error
> > output seems to round). I was therefore wondering if there is something
> I am
> > missing when it comes to deleting sites and if there is an easy way to
> > delete random sites from a 3d sample.
> >
> >
> > Thank You,
> >
> > Mitchell Greenberg
> >
> >
> > 
> > From: Joseph Weston <joseph.westo...@gmail.com>
> > Sent: Wednesday, August 30, 2017 8:33:01 PM
> > To: Mitchell Greenberg; kwant-discuss@kwant-project.org
> > Subj

Re: [Kwant] How to plot the local dos for both the scattering region and leads

[Anton Akhmerov]

Dear kuangyia,

The easiest way to resolve your problem is to add several unit cells
from the lead to the scattering region; see "add_cells" argument of
Builder.attach_lead
(https://kwant-project.org/doc/1/reference/generated/kwant.builder.Builder#kwant.builder.Builder.attach_lead)

Best,
Anton

On Mon, Sep 18, 2017 at 3:09 PM, kuangyia lee  wrote:
> Dear all,
>
> I am learning Kwant through tutorial. My toy model is a two-terminal quantum
> wire system based on the square lattice. I want to plot the spatial local
> dos for both the scattering region and leads. However, it does not work for
> leads part. The following is the code that presents my problem.  I am
> wondering whether it is possible to solve this problem with Kwant.
>
> Thanks a lot for any help.
>
> Kind regards,
>
> kuangyia
>
> =
> #!/usr/bin/env python3
>
> import kwant
> import numpy as np
> import matplotlib.pyplot as plt
>
> def make_system(a=1, t=1.0, W=20, L=50, L_well=10):
> # Start with an empty tight-binding system and a single square lattice.
> # `a` is the lattice constant (by default set to 1 for simplicity.
> lat = kwant.lattice.square(a)
> sys = kwant.Builder()
>
>  Define the scattering region and its potential profiel. 
> def potential(site, pot):
> (x, y) = site.pos
> if (L - L_well) / 2 < x < (L + L_well) / 2:
> return pot
> else:
> return 0
>
> def onsite(site, pot=0):
> return 4 * t + potential(site, pot)
>
> sys[(lat(x, y) for x in range(L) for y in range(W))] = onsite
> sys[lat.neighbors()] = -t
>
>  Define and attach the leads. 
> lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
> lead[(lat(0, j) for j in range(W))] = 4 * t
> lead[lat.neighbors()] = -t
> sys.attach_lead(lead)
> sys.attach_lead(lead.reversed())
>
> return sys
>
> def main():
>
> sys = make_system()
>
> # Check that the system looks as intended.
> kwant.plot(sys)
>
> # Finalize the system.
> sys = sys.finalized()
>
> # Calculate ldos at a given energy
> well_depth = 0.55
> ldos = kwant.ldos(sys, energy=0.25, args=[-well_depth])
> kwant.plotter.map(sys, ldos, num_lead_cells=10)
>
> if __name__ == '__main__':
> main()
>

Re: [Kwant] [kwant] lattice constant and the density of states in four-band Hamiltonian

[Anton Akhmerov]

Dear Gongwei,

Please take a look at the discretizer module introduced in Kwant 1.3 for
dealing with tight binging approximations of continuum Hamiltonians. It
should save you a lot of work and debugging.

As far as the reproducing the results of the reference, I fear it is beyond
the scope of this mailing list, and I recommend you to approach the authors
directly.

Best,
Anton

On Fri, Sep 8, 2017, 07:35 Gongwei_Hu  wrote:

> Hello everyone!
>
>
>
> I want to solve a quantum transport question about the quantum spin Hall
> bar. The Hamiltonian of my system is
>
> The tight-binding model can be written as follow ( is the lattice
> constant):
>
> Onsite item:;   Hopping-x item: ;   Hopping-y item:
>
> The matrix:
>
> My Question:
>
> (1) The first one is that I am not sure whether the tight-binding
> model presented above is right for the case of , such as , especially for
> the hopping –x and –y item.
>
> (2) In order to calculate the density distribution of electron and
> the conductance, the code presented below cannot exactly match the results
> in the paper [
> https://journals.aps.org/prb/abstract/10.1103/PhysRevB.83.081402]. For
> example the density of states calculated by Kwant near the edge is much
> thick compared with the paper results. The conductance also cannot give the
> same result, such as the conductance is 1.1e^2/h for my result but 2.0e^2/h
> for the paper at Fermi energy 10meV.
>
>
>
> My kwant code is
>
> import math
>
>
>
> import numpy
>
> from numpy import *
>
> from matplotlib import pyplot
>
>
>
> import kwant
>
> import tinyarray
>
>
>
> sigma_0 = tinyarray.array([[1, 0, 0, 0], [0, 1, 0 ,0], [0, 0, 1, 0], [0,
> 0, 0, 1]])
>
> sigma_1 = tinyarray.array([[0, -1, 0, 0], [-1, 0, 0 ,0], [0, 0, 0, 1], [0,
> 0, 1, 0]])
>
> sigma_2 = tinyarray.array([[1, 0, 0, 0], [0, -1, 0 ,0], [0, 0, 1, 0], [0,
> 0, 0, -1]])
>
> sigma_3 = tinyarray.array([[0, -1, 0, 0], [1, 0, 0 ,0], [0, 0, 0, -1], [0,
> 0, 1, 0]])
>
>
>
>
>
> def make_system(W_QPC):
>
> L1=70
>
> W1=60
>
> L2=70
>
> W=150
>
> A=364.5
>
> B=-686
>
> C=0
>
> D=-512
>
> M=-10
>
> a=2
>
> b=1
>
> def central_region(pos):
>
> x, y = pos
>
> return -(L1 +W1 + L2)/2 < x < (L1 +W1 + L2)/2 and \
>
> abs(y) < (x <= L1-(L1 +W1 + L2)/2) * abs(0.25 * (2 * W-W_QPC)
> * (math.sin((x+(L1 +W1 + L2)/2)*math.pi/L1 - math.pi/2) + 1)-W)+ (x >
> L1-(L1 +W1 + L2)/2 and x <= (L1+W1-(L1 +W1 + L2)/2)) * W_QPC/2 + (x > (L1 +
> W1-(L1 +W1 + L2)/2) and x < (L1+W1+L2-(L1 +W1 + L2)/2)) * abs(0.25 * (2 *
> W-W_QPC) * (math.sin((x-W1/2)*math.pi/L1 + math.pi/2) + 1)-W)
>
>
>
> lat = kwant.lattice.square(a)
>
> sys = kwant.Builder()
>
>
>
> sys[lat.shape(central_region, (0, 0))] = (C-4*D/a**2)*sigma_0 +
> (M-4*B/a**2)*sigma_2
>
> sys[kwant.builder.HoppingKind((b, 0), lat, lat)] = D * sigma_0/a**2+B
> * sigma_2/a**2+1j * A * sigma_1/a
>
> sys[kwant.builder.HoppingKind((0, b), lat, lat)] = D * sigma_0/a**2+B
> * sigma_2/a**2+ A * sigma_3/a
>
>
>
> sym = kwant.TranslationalSymmetry((-a, 0))
>
> lead = kwant.Builder(sym)
>
> lead[(lat(0, y) for y in range(-int(W/2) + 1, int(W/2)))] =
> (C-4*D/a**2)*sigma_0 + (M-4*B/a**2)*sigma_2
>
> lead[kwant.builder.HoppingKind((b, 0), lat, lat)] = D * sigma_0/a**2+B
> * sigma_2/a**2+1j * A * sigma_1/a
>
> lead[kwant.builder.HoppingKind((0, b), lat, lat)] = D * sigma_0/a**2+B
> * sigma_2/a**2+ A * sigma_3/a
>
>
>
> sys.attach_lead(lead)
>
> sys.attach_lead(lead.reversed())
>
>
>
> return sys.finalized()
>
>
>
> def density(sys, energy):
>
> psi = kwant.wave_function(sys, energy)
>
> for i in range(2):
>
>psi_5 = abs(psi(0)[i]**2)#.sum(axis=0)
>
> #psi_5 = abs(psi(0)[5]**2)
>
> #numpy.savetxt("psi_n.txt",psi_n);
>
>A1, A2 = psi_5[::4], psi_5[1::4]#, psi_n[2::4], psi_n[3::4]
>
>D = A1+A2
>
>numpy.savetxt("D_now.txt",D);
>
>kwant.plotter.map(sys, D)
>
>
>
> def plot_conductance(WW_QPC, energy):
>
> # Compute conductance
>
> data = []
>
> for w_QPC in WW_QPC:
>
> sys = make_system(W_QPC=w_QPC)
>
> #sys = sys.finalized()
>
> smatrix = kwant.smatrix(sys, energy)
>
> data.append(smatrix.transmission(1, 0))
>
> numpy.savetxt("WW_QPC.txt",WW_QPC)
>
> numpy.savetxt("data.txt",data)
>
> pyplot.figure()
>
> pyplot.plot(WW_QPC, data)
>
> pyplot.xlabel("W_QPC [nm]")
>
> pyplot.ylabel("conductance [e^2/h]")
>
> pyplot.show()
>
>
>
> sys = make_system(W_QPC=50)
>
> kwant.plot(sys)
>
> d = density(sys, energy=10)
>
> plot_conductance(WW_QPC = [50*i+50 for i in range(4)], energy = 10)
>
>
>
> #numpy.savetxt("D.txt",d);
>
> #kwant.plotter.map(sys, d)
>
>
>
>
>
>
>

Re: [Kwant] up/down spin conductance from Greens function like smatrix

[Anton Akhmerov]

Dear Sudin,

It is definitely possible to compute all the desired properties using the
Green's function formalism. Moreover since the Green's function is in the
basis that you specify when providing the Hamiltonian, you don't even need
any special machinery, and you can just apply the definition on your own.

Regarding the noise: there is indeed only a reference implementation of a
2-terminal shot noise. We would welcome contributions from the community
generalising that code to multiterminal situation (I believe the spin shot
noise is a particular case of multiterminal shot noise).

Best,
Anton

Best,
Anton

On Fri, Sep 15, 2017, 09:12 SUDIN GANGULY  wrote:

> Dear Sir,
>
>
> In the new version of Kwant, it is possible to calculate the up and down
> spin  conductance. In order to calculate the different components of the
> spin conductance (G_x,G_y,G_z), Greens function can be used as shown in
> the following mail-archive 
> link:https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00483.html
>
> If I want to calculate the quantities like G_x^up or G_x^down, is it
> possible to calculate these with Greens function like the smatrix.
>
> Another question I want to ask (forgive me if it is out of the context).
> We can calculate the zero temperature charge shot-noise using the
> kwant.physics.two_terminal_shotnoise. Is there any option to calculate the
> spin shot-noise?
>
> With Regards,
> Sudin Ganguly
>
>

Re: [Kwant] Understanding the scaling of the energy and wavenumber units of measurement

[Anton Akhmerov]

vector, but I do it in units of eV and nm^{-1}...
>
>
> The dirac cone is correct, as I have compared even with published papers,
> but the dispersion for the armchair is incorrect (should be contained
> within the dirac cone for such low energies).
>
> I am sorry that my issue is not solved after your answer, but I would
> appreciate any input that can make me understand why the armchair
> dispersion is not contained within the Dirac cone...
>
> Best,
> George
>
> On 10/09/2017 03:51, Anton Akhmerov wrote:
>
> Dear George,
>
> What does not make any assumptions about the units, leaving the
> interpretation entirely up to you.
>
>
> Which is the unit of measurement of length? I always thought that it is 
> simply `1`! If I write `honeycomb(32.123)` is the unit of measurement of 
> length still `1` or does it become 32.123, irrespectively of how many 
> nanometers 1 actually corresponds to? When I get the lattice sites, I know 
> that they are assumed (i,j) TIMES the lattice constant, so the actual unit of 
> distance  is still (i,j)*32.123*1, and `1` corresponds to whatever units I 
> choose.
>
> The length argument to the honeycomb lattice is the lattice constant.
> Creating a lattice does not change the units of length. The site
> positions are then calculated according to the lattice definition, so
> a site (i, j) belonging square lattice with a lattice constant a would
> have coordinates (a*i, a*j).
>
>
> Now, is the wavenumber represented in units of `1`, or in units of 1/32.123? 
> If I have a value of `k=3` and my `1` corresponds to 1 nanometer, does this 
> value of `k` correspond to `3`nm^{-1} or to `3/32.123` nm^{-1} ?
>
> kwant.physics.Bands and kwant.plotter.Bands calculate the spectrum of
> a system with a 1D translation invariance as a function of its lattice
> wave vector (so, these band structures always have a period of 2*pi).
> In other words, they measure momentum in units of inverse period of
> the system whose band structure you are calculating. This is different
> from kwant.wraparound.plot_2d_bands (introduced in kwant v1.3) that
> uses units of inverse length.
>
>
> What is the input given to the function `bands(k)`, obtained from 
> `kwant.physics.Bands` ? Does it have to always be within the interval -pi/2 
> to pi/2 , meaning that it is always expected in units of 1/lattice const. 
> irrespectively of which are the units of measurement of the wavenumber?
>
> The input to bands does not have to be within that interval, but it is
> measured in the units of the inverse period of the system, indeed.
>
>
> What is the unit of energy, that is returned by the function `bands()`? Does 
> it depend on the value that I give to the hopping, or it depends on how I 
> define `1` ? E.g. when I write `armchair[glat.neighbors()] = -1/sf` is the 
> unit of energy internally defined to be `1` or `1/sf` ??? It cannot possibly 
> be the latter, because I could also assign some number to the on-site 
> potential. So it must be the first, right? But then, how do you answer my 
> question number 4:
>
> Given that I know the length of the lattice constant in nanometers,
> and I know the energy of the hopping value in eV, how can I get the
> bands in eV versus nm^{-1} ?
>
> The unit of energy everywhere is the same as the one you use when
> defining the Hamiltonian.
>
> I hope the above answers help.
>
> Best,
> Anton
>
> On Fri, Sep 8, 2017 at 3:12 PM, George Datseris<george.datse...@ds.mpg.de> 
> <george.datse...@ds.mpg.de> wrote:
>
> Hello,
>
> I am having significant trouble understanding which are the units of
> measurement of the energy and the wavenumber when using kwant and when using
> different lattice constants.
>
> I am trying to reproduce the results of the (probably well-known paper)
> "Scalable Tight Binding Model for Graphene", Phys. Rev. Lett. 114, 036601
> (2015).
>
> I have read the answers to another question
> (https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00069.html)
> but I simply cannot understand which is the actual unit of measurement used
> internally so I can translate those to nanometers and electronVolts. Even
> though I have realized how I can translate internal wavenumbers into
> nanometers^{-1}, I cannot understand how and why to use properly
> kwants.physics.Bands().
>
> Here is the minimal code that captures my problems:
> ```python
> import kwant
> from matplotlib import pyplot
>
> a = 0.142 #inter-carbon distance in nm
> lc = sqrt(3)*a #lattice constant in nm
> t = 2.8 #hopping in eV
>
> W = 200 #width in UNSCALED lattice constants, as if sf=1. (multiply by lc to
> get it in nm)
>
> def bandstructur

Re: [Kwant] Understanding the scaling of the energy and wavenumber units of measurement

[Anton Akhmerov]

Dear George,

What does not make any assumptions about the units, leaving the
interpretation entirely up to you.

> Which is the unit of measurement of length? I always thought that it is 
> simply `1`! If I write `honeycomb(32.123)` is the unit of measurement of 
> length still `1` or does it become 32.123, irrespectively of how many 
> nanometers 1 actually corresponds to? When I get the lattice sites, I know 
> that they are assumed (i,j) TIMES the lattice constant, so the actual unit of 
> distance  is still (i,j)*32.123*1, and `1` corresponds to whatever units I 
> choose.

The length argument to the honeycomb lattice is the lattice constant.
Creating a lattice does not change the units of length. The site
positions are then calculated according to the lattice definition, so
a site (i, j) belonging square lattice with a lattice constant a would
have coordinates (a*i, a*j).

> Now, is the wavenumber represented in units of `1`, or in units of 1/32.123? 
> If I have a value of `k=3` and my `1` corresponds to 1 nanometer, does this 
> value of `k` correspond to `3`nm^{-1} or to `3/32.123` nm^{-1} ?

kwant.physics.Bands and kwant.plotter.Bands calculate the spectrum of
a system with a 1D translation invariance as a function of its lattice
wave vector (so, these band structures always have a period of 2*pi).
In other words, they measure momentum in units of inverse period of
the system whose band structure you are calculating. This is different
from kwant.wraparound.plot_2d_bands (introduced in kwant v1.3) that
uses units of inverse length.

> What is the input given to the function `bands(k)`, obtained from 
> `kwant.physics.Bands` ? Does it have to always be within the interval -pi/2 
> to pi/2 , meaning that it is always expected in units of 1/lattice const. 
> irrespectively of which are the units of measurement of the wavenumber?

The input to bands does not have to be within that interval, but it is
measured in the units of the inverse period of the system, indeed.

> What is the unit of energy, that is returned by the function `bands()`? Does 
> it depend on the value that I give to the hopping, or it depends on how I 
> define `1` ? E.g. when I write `armchair[glat.neighbors()] = -1/sf` is the 
> unit of energy internally defined to be `1` or `1/sf` ??? It cannot possibly 
> be the latter, because I could also assign some number to the on-site 
> potential. So it must be the first, right? But then, how do you answer my 
> question number 4:
Given that I know the length of the lattice constant in nanometers,
and I know the energy of the hopping value in eV, how can I get the
bands in eV versus nm^{-1} ?

The unit of energy everywhere is the same as the one you use when
defining the Hamiltonian.

I hope the above answers help.

Best,
Anton

On Fri, Sep 8, 2017 at 3:12 PM, George Datseris
 wrote:
> Hello,
>
> I am having significant trouble understanding which are the units of
> measurement of the energy and the wavenumber when using kwant and when using
> different lattice constants.
>
> I am trying to reproduce the results of the (probably well-known paper)
> "Scalable Tight Binding Model for Graphene", Phys. Rev. Lett. 114, 036601
> (2015).
>
> I have read the answers to another question
> (https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00069.html)
> but I simply cannot understand which is the actual unit of measurement used
> internally so I can translate those to nanometers and electronVolts. Even
> though I have realized how I can translate internal wavenumbers into
> nanometers^{-1}, I cannot understand how and why to use properly
> kwants.physics.Bands().
>
> Here is the minimal code that captures my problems:
> ```python
> import kwant
> from matplotlib import pyplot
>
> a = 0.142 #inter-carbon distance in nm
> lc = sqrt(3)*a #lattice constant in nm
> t = 2.8 #hopping in eV
>
> W = 200 #width in UNSCALED lattice constants, as if sf=1. (multiply by lc to
> get it in nm)
>
> def bandstructure(sf):  #sf = scaling factor
>
> glat = kwant.lattice.honeycomb(sf) #sf*la for lattice constant in nm
> A, B = glat.sublattices
> # create an armchair lead
> def leadshape(pos): #set leads width
> x, y = pos
> return (0 <= x < W)
>
> sym_lead_vertical = kwant.TranslationalSymmetry(glat.vec((-1*sf,2*sf)))
> armchair = kwant.Builder(sym_lead_vertical)
> armchair[glat.shape(leadshape, (0, 0))] = 0  #onsite energy is 0 for
> Dirac
> armchair[glat.neighbors()] = -1/sf   #hopping. use t/sf for value in eV
>
> #function that gives the bands (energies) at a given wave vector:
> bands = kwant.physics.Bands(armchair.finalized())
>
> #The wavenumbers are measured in units of 1/lattice_const.
> #that is why they only need to go from -pi to pi (Bloch):
> wavenums = np.linspace(-pi/10 ,pi/10, 41)
>
> #I now want to get the energies, measured in eV. How
> energies = [bands(k/sf)*(t/sf) for k 

Re: [Kwant] superconductivity and current conservation

[Anton Akhmerov]

Hi Tibor,

If you want to simulate a floating superconductor you only need to do an
extra step.
1. Fix the voltages of all normal leads.
2. Assuming a certain voltage of a superconductor, compute current that is
drained by it.
3. Numerically find a voltage when this current is zero by using any
root-finding algorithm.

Conceptually it isn't too hard, but it is somewhat more costly numerically.

Also of course, if you have more than one superconductor and they are at
different voltages, then this problem becomes much harder since there is no
quasiparticle energy conservation anymore.

Best,
Anton

On Thu, Aug 31, 2017 at 7:22 PM, Tibor Sekera <tibor.sek...@unibas.ch>
wrote:

> Hi Anton,
>
> I understand. Is this a feature of KWANT (and the way superconductivity is
> implemented) or of the formula for conductance?
>
> Then, for example, with KWANT we cannot calculate properly the conductance
> of a NSN system, where S is finite (not grounded) superconductor
> and N are leads. Because if we make S region very large we can get G_00
> finite while G_10=0.
>
> In
> https://journals.aps.org/prb/abstract/10.1103/PhysRevB.53.16390
> they claim that Eq.(31) is valid also for a 'floating' superconductor (not
> grounded finite piece), like in their Fig.2(b).
>
> In other words, how to treat properly with KWANT such a case, when we want
> to have finite piece of superconductor, that is not grounded?
>
> Best,
> Tibor
> --
> *From:* anton.akhme...@gmail.com [anton.akhme...@gmail.com] on behalf of
> Anton Akhmerov [anton.akhmerov...@gmail.com]
> *Sent:* Thursday, August 31, 2017 6:58 PM
> *To:* Tibor Sekera
> *Cc:* kwant-discuss@kwant-project.org
> *Subject:* Re: [Kwant] superconductivity and current conservation
>
> Hi Tibor,
>
> As soon as there is a superconducting region of any size, the charge can
> leave away as supercurrent, while the current carried by the quasiparticles
> isn't conserved anymore. Essentially all the superconductors in the system,
> both finite and infinite are an extra lead.
>
> I hope this answers your question,
> Anton
>
> On Thu, Aug 31, 2017 at 6:25 PM, Tibor Sekera <tibor.sek...@unibas.ch>
> wrote:
>
>> Dear all,
>>
>> in the tutorial
>> https://kwant-project.org/doc/1/tutorial/superconductors
>> one calculates the conductance from left lead to right lead, G_10, as
>> G_00 = e^2/h (N-R_ee+R_he)
>> assuming that, due to the current conservation,
>> G_10 = G_00.
>>
>> My question is, what happens if we detach the superconducting lead L1.
>> Intuitively, I would expect
>> G_00 = 0. But using the formula above one obtains non-zero values, and
>> making the finite superconducting
>> region (the part of the scattering region) large enough, the subgap
>> transport is exactly as if
>> the lead L1 was attached. This means that the current conservation is
>> violated, and something is wrong.
>>
>> Below is slightly modified code from the tutorial to generate a plot
>> comparing G_00 for a system with and
>> without lead L1. You can see, that if parameter L (the superconducting
>> part of the scattering region) is large
>> enough, the subgap transport is the same whether L1 is attached or not.
>>
>> Why is this happening, can someone point out the error?
>>
>> Best,
>> Tibor
>>
>>
>>
>> import kwant
>> import tinyarray
>> import numpy as np
>> %matplotlib inline
>> from matplotlib import pyplot as plt
>>
>> tau_x = tinyarray.array([[0, 1], [1, 0]])
>> tau_y = tinyarray.array([[0, -1j], [1j, 0]])
>> tau_z = tinyarray.array([[1, 0], [0, -1]])
>>
>> L=80
>> barrier=1.5
>> W=10
>>
>> def make_system(a=1, W=W, L=L, barrier=barrier, barrierpos=(3, 4),
>> mu=0.4, Delta=0.1, Deltapos=4, t=1.0, phs=True):
>> lat = kwant.lattice.square(norbs=2)
>> syst = kwant.Builder()
>>
>> syst[(lat(x, y) for x in range(0,Deltapos) for y in range(W))] = (4 *
>> t - mu) * tau_z
>> syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4
>> * t - mu) * tau_z + Delta * tau_x
>>
>> # The tunnel barrier
>> syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1])
>>  for y in range(W))] = (4 * t + barrier - mu) * tau_z
>>
>> # Hoppings
>> syst[lat.neighbors()] = -t * tau_z
>> sym_left = kwant.TranslationalSymmetry((-a, 0))
>> lead0 = kwant.Builder(sym_left, conservation_law=-tau_z,
>> particle_hole=tau_y)
>> lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z
>> lead0[lat.neighb

Re: [Kwant] superconductivity and current conservation

[Anton Akhmerov]

Hi Tibor,

As soon as there is a superconducting region of any size, the charge can
leave away as supercurrent, while the current carried by the quasiparticles
isn't conserved anymore. Essentially all the superconductors in the system,
both finite and infinite are an extra lead.

I hope this answers your question,
Anton

On Thu, Aug 31, 2017 at 6:25 PM, Tibor Sekera 
wrote:

> Dear all,
>
> in the tutorial
> https://kwant-project.org/doc/1/tutorial/superconductors
> one calculates the conductance from left lead to right lead, G_10, as
> G_00 = e^2/h (N-R_ee+R_he)
> assuming that, due to the current conservation,
> G_10 = G_00.
>
> My question is, what happens if we detach the superconducting lead L1.
> Intuitively, I would expect
> G_00 = 0. But using the formula above one obtains non-zero values, and
> making the finite superconducting
> region (the part of the scattering region) large enough, the subgap
> transport is exactly as if
> the lead L1 was attached. This means that the current conservation is
> violated, and something is wrong.
>
> Below is slightly modified code from the tutorial to generate a plot
> comparing G_00 for a system with and
> without lead L1. You can see, that if parameter L (the superconducting
> part of the scattering region) is large
> enough, the subgap transport is the same whether L1 is attached or not.
>
> Why is this happening, can someone point out the error?
>
> Best,
> Tibor
>
>
>
> import kwant
> import tinyarray
> import numpy as np
> %matplotlib inline
> from matplotlib import pyplot as plt
>
> tau_x = tinyarray.array([[0, 1], [1, 0]])
> tau_y = tinyarray.array([[0, -1j], [1j, 0]])
> tau_z = tinyarray.array([[1, 0], [0, -1]])
>
> L=80
> barrier=1.5
> W=10
>
> def make_system(a=1, W=W, L=L, barrier=barrier, barrierpos=(3, 4),
> mu=0.4, Delta=0.1, Deltapos=4, t=1.0, phs=True):
> lat = kwant.lattice.square(norbs=2)
> syst = kwant.Builder()
>
> syst[(lat(x, y) for x in range(0,Deltapos) for y in range(W))] = (4 *
> t - mu) * tau_z
> syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 *
> t - mu) * tau_z + Delta * tau_x
>
> # The tunnel barrier
> syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1])
>  for y in range(W))] = (4 * t + barrier - mu) * tau_z
>
> # Hoppings
> syst[lat.neighbors()] = -t * tau_z
> sym_left = kwant.TranslationalSymmetry((-a, 0))
> lead0 = kwant.Builder(sym_left, conservation_law=-tau_z,
> particle_hole=tau_y)
> lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z
> lead0[lat.neighbors()] = -t * tau_z
> sym_right = kwant.TranslationalSymmetry((a, 0))
> lead1 = kwant.Builder(sym_right)
> lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta *
> tau_x
> lead1[lat.neighbors()] = -t * tau_z
> syst.attach_lead(lead0)
> syst.attach_lead(lead1)
> return syst
>
> syst = make_system()
> kwant.plot(syst)
> syst = syst.finalized()
>
> data = []
> energies=[0.002 * i for i in range(0,100)]
> for energy in energies:
> smatrix = kwant.smatrix(syst, energy)
> # Conductance is N - R_ee + R_he
> data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] -
> smatrix.transmission((0, 0), (0, 0)) +
> smatrix.transmission((0, 1), (0, 0)))
>
> syst = make_system()
> del syst.leads[-1]
> # Check that the system looks as intended.
> kwant.plot(syst)
> # Finalize the system.
> syst = syst.finalized()
> data2 = []
> for energy in energies:
> smatrix = kwant.smatrix(syst, energy)
> # Conductance is N - R_ee + R_he
> data2.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] -
> smatrix.transmission((0, 0), (0, 0)) +
> smatrix.transmission((0, 1), (0, 0)))
>
> fig , (ax0) = plt.subplots(1, 1, sharey = True, figsize = [5, 5])
> ax0.plot(energies, data, label = 'with lead 1')
> ax0.plot(energies, data2, label = 'without lead 1')
> ax0.set_xlabel("energy [t]")
> ax0.set_ylabel("conductance [e^2/h]")
> ax0.set_ylim(0,3.5)
> ax0.legend()
> fig.savefig('NS_conductance.png')
>

Re: [Kwant] 2D Current Density Streamplot of 3D System

[Anton Akhmerov]

> Yep okay that makes sense now, sorry that didn't click. One final question I
> would have then is that the field vector array has a z index of 28, whereas
> my system has a height of 52 with 10 layers in it, is this due to the
> prefactor that I might need to redefine?

No, the resolution is controlled by other arguments of interpolate
current, please consult the documentation.

Best,
Anton

> Thank You,
>
> Mitchell Greenberg
>
> 
> From: anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of Anton
> Akhmerov <anton.akhmerov...@gmail.com>
> Sent: Sunday, August 27, 2017 7:30:44 PM
>
> To: Mitchell Greenberg
> Cc: kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>
> Hi MItchel,
>
>> Sorry this stuff is all a bit beyond me. If I apply current[:, :, 12, :2].
>> am I essentially taking a slice at z = 0? If so how would I modify that
>> expression to look at a layer where z = 20?
>
> That would be z = 12 (so you'd need to change the third index). If you
> want to understand numpy slicing better, take a look here:
> https://arxiv.org/abs/1102.1523
>
> Best,
> Anton
>
>> Sorry for the trouble,
>>
>> Mitchell Greenberg
>>
>> 
>> From: anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of
>> Anton
>> Akhmerov <anton.akhmerov...@gmail.com>
>> Sent: Sunday, August 27, 2017 6:34:58 PM
>>
>> To: Mitchell Greenberg
>> Cc: kwant-discuss@kwant-project.org
>> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>>
>> Hi Mitchell,
>>
>>> The shape of the field array that gets returned is:
>>>
>>> (38, 28, 28, 3)
>>>
>>> From looking at the source code the 3 on the end seems to just be for
>>> whatever dimension your system is, so I assume it's not a simple fix as
>>> it
>>> seems kwants 'streamplot' wants a  3d array and not a 4d array.
>>
>> Since current is a vector and not a scalar, you get a vector field
>> evaluated over a homogeneous grid. And of course now it is a
>> 3D-vector, and not a 2D vector. So for example you could get rid of
>> the z-component entirely by selecting current[:, :, 12, :2].
>>
>> Another interesting possibility that I suggest to explore is to
>> somehow use ipyvolume (see ipyvolume.readthedocs.io) for rendering the
>> 3D current distribution.
>>
>> Cheers,
>> Anton
>>
>>>
>>>
>>> Any advice would be greatly appreciated.
>>>
>>>
>>> Thank You,
>>>
>>> Mitchell Greenberg
>>>
>>> 
>>> From: anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of
>>> Anton
>>> Akhmerov <anton.akhmerov...@gmail.com>
>>> Sent: Saturday, August 26, 2017 8:49:08 PM
>>>
>>> To: Mitchell Greenberg
>>> Cc: kwant-discuss@kwant-project.org
>>> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>>>
>>> Hi Mitchell,
>>>
>>> The only problem that occurs if you remove the check is the
>>> normalization of the current density; in 3D you'd need to redo the
>>> integration of the smoothing function to get the right prefactor. If
>>> that doesn't bother you, just remove the check (or you could compute
>>> the correct prefactor).
>>>
>>> Best,
>>> Anton
>>>
>>> On Sat, Aug 26, 2017 at 3:30 AM, Mitchell Greenberg
>>> <mitchell.greenb...@my.jcu.edu.au> wrote:
>>>> Hi Anton,
>>>>
>>>>
>>>> Thanks for the quick reply, when I try to call
>>>> kwant.plotter.interpolate_current(syst,current) for my 3D system, I
>>>> encounter the following error:
>>>>
>>>>
>>>> 'interpolate_current' only works for 2D systems.
>>>>
>>>>
>>>> The error triggers at the line:
>>>>
>>>>
>>>> # TODO: Generalize 'factor' prefactor to arbitrary dimensions and
>>>> remove
>>>> #   this check. This check is done here to keep changes local
>>>> if dim != 2:
>>>> raise ValueError("'interpolate_current' only works for 2D
>>>> systems.")
>>>>
>>>>
>>>> I went and had a look at the source code and it seems that this check is
>>>> here 'to keep changes local' and that a 'TODO' is to generaliz

Re: [Kwant] 2D Current Density Streamplot of 3D System

[Anton Akhmerov]

Hi MItchel,

> Sorry this stuff is all a bit beyond me. If I apply current[:, :, 12, :2].
> am I essentially taking a slice at z = 0? If so how would I modify that
> expression to look at a layer where z = 20?

That would be z = 12 (so you'd need to change the third index). If you
want to understand numpy slicing better, take a look here:
https://arxiv.org/abs/1102.1523

Best,
Anton

> Sorry for the trouble,
>
> Mitchell Greenberg
>
> 
> From: anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of Anton
> Akhmerov <anton.akhmerov...@gmail.com>
> Sent: Sunday, August 27, 2017 6:34:58 PM
>
> To: Mitchell Greenberg
> Cc: kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>
> Hi Mitchell,
>
>> The shape of the field array that gets returned is:
>>
>> (38, 28, 28, 3)
>>
>> From looking at the source code the 3 on the end seems to just be for
>> whatever dimension your system is, so I assume it's not a simple fix as it
>> seems kwants 'streamplot' wants a  3d array and not a 4d array.
>
> Since current is a vector and not a scalar, you get a vector field
> evaluated over a homogeneous grid. And of course now it is a
> 3D-vector, and not a 2D vector. So for example you could get rid of
> the z-component entirely by selecting current[:, :, 12, :2].
>
> Another interesting possibility that I suggest to explore is to
> somehow use ipyvolume (see ipyvolume.readthedocs.io) for rendering the
> 3D current distribution.
>
> Cheers,
> Anton
>
>>
>>
>> Any advice would be greatly appreciated.
>>
>>
>> Thank You,
>>
>> Mitchell Greenberg
>>
>> 
>> From: anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of
>> Anton
>> Akhmerov <anton.akhmerov...@gmail.com>
>> Sent: Saturday, August 26, 2017 8:49:08 PM
>>
>> To: Mitchell Greenberg
>> Cc: kwant-discuss@kwant-project.org
>> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>>
>> Hi Mitchell,
>>
>> The only problem that occurs if you remove the check is the
>> normalization of the current density; in 3D you'd need to redo the
>> integration of the smoothing function to get the right prefactor. If
>> that doesn't bother you, just remove the check (or you could compute
>> the correct prefactor).
>>
>> Best,
>> Anton
>>
>> On Sat, Aug 26, 2017 at 3:30 AM, Mitchell Greenberg
>> <mitchell.greenb...@my.jcu.edu.au> wrote:
>>> Hi Anton,
>>>
>>>
>>> Thanks for the quick reply, when I try to call
>>> kwant.plotter.interpolate_current(syst,current) for my 3D system, I
>>> encounter the following error:
>>>
>>>
>>> 'interpolate_current' only works for 2D systems.
>>>
>>>
>>> The error triggers at the line:
>>>
>>>
>>> # TODO: Generalize 'factor' prefactor to arbitrary dimensions and
>>> remove
>>> #   this check. This check is done here to keep changes local
>>> if dim != 2:
>>> raise ValueError("'interpolate_current' only works for 2D
>>> systems.")
>>>
>>>
>>> I went and had a look at the source code and it seems that this check is
>>> here 'to keep changes local' and that a 'TODO' is to generalize the
>>> prefactor to arbitrary dimensions, can I just   remove this check or will
>>> that cause issues?
>>>
>>>
>>> If this cannot be done I was looking at current density calculations on
>>> the
>>> forums last night and I saw there was code for a manual calculation of
>>> current density, I also saw on the forums that we can get the 'slice' of
>>> the
>>> wavefunctions through:
>>>
>>>
>>> wavefunction = wf(0)
>>>
>>> def in_sheet(site):
>>> return site.pos[1] == 0  # site in x-z plane, y position == 0
>>>
>>>
>>> For my code I'd switch to site.pos[2] == 0 since we want z axis slices
>>>
>>>
>>> sheet = [idx for idx, site in enumerate(sysf.sites) if
>>> in_sheet(site)]
>>>
>>>
>>> wavefunction_in_plane = wavefunction[0][sheet]
>>>
>>>
>>> Is there a way I can use the manual calculation:
>>>
>>>
>>> I_ij = 2 핀m[ Ψ⁺_i H_ij Ψ_j ]
>>>
>>>
>>> By only evaluating the 2D section of the wavefunction and hamiltonian
&g

Re: [Kwant] 2D Current Density Streamplot of 3D System

[Anton Akhmerov]

Hi Mitchell,

> The shape of the field array that gets returned is:
>
> (38, 28, 28, 3)
>
> From looking at the source code the 3 on the end seems to just be for
> whatever dimension your system is, so I assume it's not a simple fix as it
> seems kwants 'streamplot' wants a  3d array and not a 4d array.

Since current is a vector and not a scalar, you get a vector field
evaluated over a homogeneous grid. And of course now it is a
3D-vector, and not a 2D vector. So for example you could get rid of
the z-component entirely by selecting current[:, :, 12, :2].

Another interesting possibility that I suggest to explore is to
somehow use ipyvolume (see ipyvolume.readthedocs.io) for rendering the
3D current distribution.

Cheers,
Anton

>
>
> Any advice would be greatly appreciated.
>
>
> Thank You,
>
> Mitchell Greenberg
>
> 
> From: anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of Anton
> Akhmerov <anton.akhmerov...@gmail.com>
> Sent: Saturday, August 26, 2017 8:49:08 PM
>
> To: Mitchell Greenberg
> Cc: kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>
> Hi Mitchell,
>
> The only problem that occurs if you remove the check is the
> normalization of the current density; in 3D you'd need to redo the
> integration of the smoothing function to get the right prefactor. If
> that doesn't bother you, just remove the check (or you could compute
> the correct prefactor).
>
> Best,
> Anton
>
> On Sat, Aug 26, 2017 at 3:30 AM, Mitchell Greenberg
> <mitchell.greenb...@my.jcu.edu.au> wrote:
>> Hi Anton,
>>
>>
>> Thanks for the quick reply, when I try to call
>> kwant.plotter.interpolate_current(syst,current) for my 3D system, I
>> encounter the following error:
>>
>>
>> 'interpolate_current' only works for 2D systems.
>>
>>
>> The error triggers at the line:
>>
>>
>> # TODO: Generalize 'factor' prefactor to arbitrary dimensions and
>> remove
>> #   this check. This check is done here to keep changes local
>> if dim != 2:
>> raise ValueError("'interpolate_current' only works for 2D
>> systems.")
>>
>>
>> I went and had a look at the source code and it seems that this check is
>> here 'to keep changes local' and that a 'TODO' is to generalize the
>> prefactor to arbitrary dimensions, can I just   remove this check or will
>> that cause issues?
>>
>>
>> If this cannot be done I was looking at current density calculations on
>> the
>> forums last night and I saw there was code for a manual calculation of
>> current density, I also saw on the forums that we can get the 'slice' of
>> the
>> wavefunctions through:
>>
>>
>> wavefunction = wf(0)
>>
>> def in_sheet(site):
>> return site.pos[1] == 0  # site in x-z plane, y position == 0
>>
>>
>> For my code I'd switch to site.pos[2] == 0 since we want z axis slices
>>
>>
>> sheet = [idx for idx, site in enumerate(sysf.sites) if in_sheet(site)]
>>
>>
>> wavefunction_in_plane = wavefunction[0][sheet]
>>
>>
>> Is there a way I can use the manual calculation:
>>
>>
>> I_ij = 2 핀m[ Ψ⁺_i H_ij Ψ_j ]
>>
>>
>> By only evaluating the 2D section of the wavefunction and hamiltonian that
>> we want?
>>
>>
>> Many Thanks,
>>
>> Mitchell Greenberg
>>
>> 
>> From: anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of
>> Anton
>> Akhmerov <anton.akhmerov...@gmail.com>
>> Sent: Saturday, August 26, 2017 2:35:56 AM
>> To: Mitchell Greenberg
>> Cc: kwant-discuss@kwant-project.org
>> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>>
>> Hi Mitchell,
>>
>> You can use the low level function that calculates the current
>> interpolation:
>>
>> https://kwant-project.org/doc/1/reference/generated/kwant.plotter.interpolate_current#kwant.plotter.interpolate_current
>>
>> I believe this function will work in 3D as well. Then you can slice
>> the result to get the 2D distribution that you need.
>>
>> Internally we do that and then essentially pass the output to the
>> matplotlib streamplot (or rather kwant.plotter.streamplot).
>>
>> Best,
>> Anton
>>
>> On Fri, Aug 25, 2017 at 12:15 PM, Mitchell Greenberg
>> <mitchell.greenb...@my.jcu.edu.au> wrote:
>>> Hello,
>>>
>>>
>>>
>>> I am currently completing my thesis using Kwant, per discussions with my
>>> supervisor we would like to plot the current density as a streamplot of
>>> our
>>> 3d system for different z-axis 'slices' of the system (our system has an
>>> exponentially decaying potential in the z direction, so we'd expect
>>> current
>>> to flow differently through each layer).
>>>
>>>
>>>
>>> I have tried tinkering with the new current functions implemented in 1.3
>>> but
>>> the plotting function can only plot 2d systems. Is there a way to obtain
>>> the
>>> current density streamplot for a 2d 'slice' of our 3d system.
>>>
>>>
>>>
>>> Thank You,
>>>
>>> Mitchell Greenberg
>>>
>>>

Re: [Kwant] 2D Current Density Streamplot of 3D System

[Anton Akhmerov]

Hi Mitchell,

The only problem that occurs if you remove the check is the
normalization of the current density; in 3D you'd need to redo the
integration of the smoothing function to get the right prefactor. If
that doesn't bother you, just remove the check (or you could compute
the correct prefactor).

Best,
Anton

On Sat, Aug 26, 2017 at 3:30 AM, Mitchell Greenberg
<mitchell.greenb...@my.jcu.edu.au> wrote:
> Hi Anton,
>
>
> Thanks for the quick reply, when I try to call
> kwant.plotter.interpolate_current(syst,current) for my 3D system, I
> encounter the following error:
>
>
> 'interpolate_current' only works for 2D systems.
>
>
> The error triggers at the line:
>
>
> # TODO: Generalize 'factor' prefactor to arbitrary dimensions and remove
> #   this check. This check is done here to keep changes local
> if dim != 2:
> raise ValueError("'interpolate_current' only works for 2D systems.")
>
>
> I went and had a look at the source code and it seems that this check is
> here 'to keep changes local' and that a 'TODO' is to generalize the
> prefactor to arbitrary dimensions, can I just   remove this check or will
> that cause issues?
>
>
> If this cannot be done I was looking at current density calculations on the
> forums last night and I saw there was code for a manual calculation of
> current density, I also saw on the forums that we can get the 'slice' of the
> wavefunctions through:
>
>
> wavefunction = wf(0)
>
> def in_sheet(site):
> return site.pos[1] == 0  # site in x-z plane, y position == 0
>
>
> For my code I'd switch to site.pos[2] == 0 since we want z axis slices
>
>
> sheet = [idx for idx, site in enumerate(sysf.sites) if in_sheet(site)]
>
>
> wavefunction_in_plane = wavefunction[0][sheet]
>
>
> Is there a way I can use the manual calculation:
>
>
> I_ij = 2 핀m[ Ψ⁺_i H_ij Ψ_j ]
>
>
> By only evaluating the 2D section of the wavefunction and hamiltonian that
> we want?
>
>
> Many Thanks,
>
> Mitchell Greenberg
>
> 
> From: anton.akhme...@gmail.com <anton.akhme...@gmail.com> on behalf of Anton
> Akhmerov <anton.akhmerov...@gmail.com>
> Sent: Saturday, August 26, 2017 2:35:56 AM
> To: Mitchell Greenberg
> Cc: kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System
>
> Hi Mitchell,
>
> You can use the low level function that calculates the current
> interpolation:
> https://kwant-project.org/doc/1/reference/generated/kwant.plotter.interpolate_current#kwant.plotter.interpolate_current
>
> I believe this function will work in 3D as well. Then you can slice
> the result to get the 2D distribution that you need.
>
> Internally we do that and then essentially pass the output to the
> matplotlib streamplot (or rather kwant.plotter.streamplot).
>
> Best,
> Anton
>
> On Fri, Aug 25, 2017 at 12:15 PM, Mitchell Greenberg
> <mitchell.greenb...@my.jcu.edu.au> wrote:
>> Hello,
>>
>>
>>
>> I am currently completing my thesis using Kwant, per discussions with my
>> supervisor we would like to plot the current density as a streamplot of
>> our
>> 3d system for different z-axis 'slices' of the system (our system has an
>> exponentially decaying potential in the z direction, so we'd expect
>> current
>> to flow differently through each layer).
>>
>>
>>
>> I have tried tinkering with the new current functions implemented in 1.3
>> but
>> the plotting function can only plot 2d systems. Is there a way to obtain
>> the
>> current density streamplot for a 2d 'slice' of our 3d system.
>>
>>
>>
>> Thank You,
>>
>> Mitchell Greenberg
>>
>>

Re: [Kwant] 2D Current Density Streamplot of 3D System

[Anton Akhmerov]

Hi Mitchell,

You can use the low level function that calculates the current
interpolation: 
https://kwant-project.org/doc/1/reference/generated/kwant.plotter.interpolate_current#kwant.plotter.interpolate_current

I believe this function will work in 3D as well. Then you can slice
the result to get the 2D distribution that you need.

Internally we do that and then essentially pass the output to the
matplotlib streamplot (or rather kwant.plotter.streamplot).

Best,
Anton

On Fri, Aug 25, 2017 at 12:15 PM, Mitchell Greenberg
 wrote:
> Hello,
>
>
>
> I am currently completing my thesis using Kwant, per discussions with my
> supervisor we would like to plot the current density as a streamplot of our
> 3d system for different z-axis 'slices' of the system (our system has an
> exponentially decaying potential in the z direction, so we'd expect current
> to flow differently through each layer).
>
>
>
> I have tried tinkering with the new current functions implemented in 1.3 but
> the plotting function can only plot 2d systems. Is there a way to obtain the
> current density streamplot for a 2d 'slice' of our 3d system.
>
>
>
> Thank You,
>
> Mitchell Greenberg
>
>

Re: [Kwant] Importing coordinates of the scattering region

[Anton Akhmerov]

Hi Patrik,

You can use site coordinates that don't belong to any lattice. For
that you need to define your own SiteFamily, see for example the
script below that achieves this goal.

---
from matplotlib import pyplot
import kwant
import tinyarray

class Amorphous(kwant.builder.SimpleSiteFamily):
def normalize_tag(self, tag):
return tinyarray.array(tag, float)

def pos(self, tag):
return tag


atoms = Amorphous()

syst = kwant.Builder()
sites = atoms(0.01, 0.05), atoms(1.01, -.5)
syst[sites[0]] = syst[sites[1]] = 2
syst[sites[1], sites[0]] = 1

kwant.plot(syst)
-

Of course it's now your responsibility to specify the Hamiltonian,
leads, attaching leads, etc.

Best,
Anton

On Thu, Jul 6, 2017 at 4:12 PM, Patrik Arvoy  wrote:
>
> Dear users and developers,
>
> I was wondering if I can import the coordinates of a scattering region
> without any particular symmetry and compute the conductance via kwant.
> If yes, how does one do that?
> Let's say the scattering region does not have any symmetry after optimizing
> the structure via MD or ... and I have a list of coordinates of all the
> sites.
> I appreciate any help.
>
> Regards
> Patrik

Re: [Kwant] Kwant 1.3 Installation

[Anton Akhmerov]

Dear Qingtian,

This is very weird. Can you please try running a script with the
following lines:

import sympy
from kwant.continuum.discretizer import discretize,
discretize_symbolic, build_discretized
from kwant.continuum._common import sympify, lambdify
from kwant.continuum._common import momentum_operators, position_operators

If I understand correctly, these should fail with a more informative error.

Best,
Anton

On Thu, Jun 8, 2017 at 10:30 AM, Qingtian Zhang
 wrote:
> Dear all,
> I have installed Kwant 1.3 on my computer according to the  instructions,
> but i got some error message when i ran the source code for tutorial 2.10. I
> can run all the other source codes for tutorials 2.7-2.9. The information
> for my pc is: Microsoft Windows 7， 64-bit. Python 3.6.
>
> The message I got:
> Traceback (most recent call last):
>   File "C:\Users\QingtianZhang\Desktop\discretize(1).py", line 196, in
> 
> main()
>   File "C:\Users\QingtianZhang\Desktop\discretize(1).py", line 183, in main
> template = kwant.continuum.discretize('k_x * A(x) * k_x')
>   File "C:\Python36\lib\site-packages\kwant\_common.py", line 61, in
> __getattr__
> raise RuntimeError(msg.format(self.name, self.dependencies))
> RuntimeError: 'kwant.continuum' is not available because one or more of the
> following dependencies are not installed: sympy
>
> I have downloaded 13 packages and installed all of them using pip:
> cycler-0.10.0-py2.py3-none-any
> kwant-1.3.1-cp36-cp36m-win_amd64
> matplotlib-2.0.2-cp36-cp36m-win_amd64
> numpy-1.13.0rc2+mkl-cp36-cp36m-win_amd64
> pyparsing-2.2.0-py2.py3-none-any
> pytest-3.1.1-py2.py3-none-any
> python_dateutil-2.6.0-py2.py3-none-any
> pytz-2017.2-py2.py3-none-any
> scipy-0.19.0-cp36-cp36m-win_amd64
> setuptools-35.0.2-py2.py3-none-any
> six-1.10.0-py2.py3-none-any
> sympy-1.0-py2.py3-none-any
> tinyarray-1.2.0-cp36-cp36m-win_amd64
>
> I have tried to install sympy again using:
> C:\Users\QingtianZhang>cd c:\Kwant13
> c:\Kwant13>pip3 install --no-deps sympy-1.0-py2.py3-none-any.whl
> Requirement already satisfied: sympy==1.0 from file:///C:/Kwant13/sympy-1.0-
> py2.
> py3-none-any.whl in c:\python36\lib\site-packages
>
> It shows that i have already installed sympy.
> How to fix this?
> Sincerely,
> Qingtian Zhang
>
>
>

Re: [Kwant] performance of mumps solver and sparse solver.

[Anton Akhmerov]

Dear KS,

On a linux workstation the ratio of runtimes is around 2.3 in favor of
MUMPS, so the problem seems to pertain to your installation or to the
windows installations in general. Can anyone using Kwant on windows
try the attached script? (I reduced the number of points so that it
doesn't take 1/2 hour and removed some parts not crucial to the
benchmark)

Best,
Anton

On Tue, May 16, 2017 at 7:18 PM, Prof. CHAN Kwok Sum
<apksc...@cityu.edu.hk> wrote:
> Dear Anton,
> To confirm my observation, I remove the non-essential part of the script to 
> test the mumps solver and default solver. The system is a zigzag graphene 
> ribbon, with 200X200 unit cells, there are 4 atoms in one unit cell so there 
> are 1.6x10^5 atoms. The OS is window 10, kwant is 1.2.2 and python is 3.4.4. 
> the time for calculating 20 smatrix is around 30 minutes for both default 
> solver and mumps. Please find attached the script and the runtime. The first 
> runtime is for mumps. The second runtime is for default. The kwant was 
> installed with the package prepared by Gohlke.
> For the sparse solver the time is longer, 40 minutes but there is no umfpack 
> installed in the window system. So I do not show the runtime as it may not be 
> accurate.
> KS
>
> On 16/5/2017, 5:51 PM, "anton.akhme...@gmail.com on behalf of Anton Akhmerov" 
> <anton.akhme...@gmail.com on behalf of anton.akhmerov...@gmail.com> wrote:
>
> Dear KS,
>
> > I have compared the performance of default, mumps solvers in windows and
> > Ubuntu machines. Their performance is very similar without significant 
> gain
> > when using mumps. Do I miss anything in the installation? For example, 
> not
> > installing the mumps. Similar performance is found for the sparse 
> solver.
> > The system I considered has 10^5 to 6X10^5 sites and machine has a quad 
> core
> > i7 cpu with 16G RAM. My installation procedures are those found in the
> > website for pre complled packages.
>
> This is contrary to our experience. Can you please provide more details:
> * How did you install Kwant in both cases?
> * What specific script did you run?
> * What were the execution times that you observed?
>
> Thanks,
> Anton
>
> > KS
> >
> >
> >
> > 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.
>
>
>
>
> 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.
import os
from time import time
from math import pi, sqrt, tanh
from cmath import exp


import kwant
import kwant.solvers.sparse as ksparse
import kwant.solvers.mumps as kmumps
import numpy as np


###setup directory
##curdir=os.getcwd()
##os.chdir(curdir+'/test')
print (ksparse.uses_umfpack)
#define function for magnetic field phase shift
def hopphase_x(sitei, sitej, flux):

xi, yi = sitei.pos
xj, yj = sitej.pos
return exp(-0.5*1j*flux*(xi-xj)*(yi+yj))
def hopphase_y(sitei, sitej, flux):
 xi, yi = sitei.pos
 xj, yj = sitej.pos
 return exp(0.5*1j* flux * (xi + xj) * (yi - yj))

#HIDDEN_END_dwhx
print ('start',time())
sys = kwant.Builder()
sys2=kwant.Builder()
#HIDDEN_END_goiq

# Here, we are only working with square lattices
#HIDDEN_BEGIN_suwo
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
#a = 1
# lat = kwant.lattice.honeycomb(a)
# this lattice is not a honeycomb lattice. it is like a rectangle
# this lattice gives correct zigzag and armchair edges
graphene= kwant.lattice.general([(sqrt(3), 0), (0, 1)],
   [(0,1/2), (1/(2*sqrt(3)), 0),(sqrt(3)/2, 0), (2/sqrt(3),1/2)])



def rect(pos):
x, y= pos
lx=L*sqrt(3);ly=W*1;
return 0<=x=i>=lpos: wtmp=-Lstub
for j in range(wtmp,W):
# On-site Hamiltonian
# if i<23 : 
#sys[a(i, j)] = 0. * t
#sys[b(i,j)]=0.*t
#sys[c(i,j)]=0.*t
#sys[d(i,j)]=0.*t
#  else :
sys[a(i, j)] = (0.+i*0.0)* t+shift

Re: [Kwant] performance of mumps solver and sparse solver.

[Anton Akhmerov]

Dear KS,

> I have compared the performance of default, mumps solvers in windows and
> Ubuntu machines. Their performance is very similar without significant gain
> when using mumps. Do I miss anything in the installation? For example, not
> installing the mumps. Similar performance is found for the sparse solver.
> The system I considered has 10^5 to 6X10^5 sites and machine has a quad core
> i7 cpu with 16G RAM. My installation procedures are those found in the
> website for pre complled packages.

This is contrary to our experience. Can you please provide more details:
* How did you install Kwant in both cases?
* What specific script did you run?
* What were the execution times that you observed?

Thanks,
Anton

> KS
>
>
>
> 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.

Re: [Kwant] Making 3D leads in a monoatomic lattice

[Anton Akhmerov]

Dear John,

I believe the error message is quite accurate here: (0, 0, a) is not a
lattice vector, and therefore your lattice may not have such a
translational symmetry. If you double the period of the lead symmetry, the
lead creation should work.

Best,
Anton

On Thu, Apr 13, 2017, 16:57 John Eaves  wrote:

> Good day
>
> I am trying to make a simple 3D geometry; a cylinder. Now, the only
> tutorial available for 3D is
> https://kwant-project.org/doc/1.0/tutorial/tutorial6#d-example-zincblende-structure
> so there is a chance my problems come from that.
>
> The thing is that I cannot seem to attach any leads. The scattering region
> code works just fine, but when I try to attach the leads I get the error
>
> Site family  [1.0 1.0 0.0], origin [0.0 0.0 0.0]> does not have commensurate periods with 
> symmetry .
> The error is not completely opague of course; something seems to be 
> problematic about the lattice vectors and the proposed symmetry. For that it 
> is perhaps good to include the full code:
>
> def make_system(a=1, t=1.0, W=10, L=5, r2=20):
> lat = kwant.lattice.general([(0, a, a), (a, 0, a), (a, a, 0)])
> sys = kwant.Builder()
>
> def ring(pos):
> (x, y,z) = pos
> rsq = x ** 2 + y ** 2
> return rsq < r2 ** 2 and 0 <= z < L
>
> sys[lat.shape(ring, (1, 1,1))] = 4 * t
> sys[lat.neighbors()] = -t
>
> sym_lead = kwant.TranslationalSymmetry((0, 0,-a))
> lead = kwant.Builder(sym_lead)
>
> def lead_shape(pos):
> (x, y,z) = pos
> rsq = x ** 2 + y ** 2
> return rsq < r2 ** 2
>
> lead[lat.shape(lead_shape, (0, 0,0))] = 4 * t
> lead[lat.neighbors()] = -t
>
> sys.attach_lead(lead)
> sys.attach_lead(lead.reversed())
>
> return sys
>
>
>
> I personally do not see what is wrong with this, especially as the scattering 
> region looks fine without the leads. According to the error there is 
> something wrong with the lattice vectors themselves? Should I change them in 
> some way?
> I don't care about any particular structure here really, I am just looking 
> for the 3D equivalent of the square lattice and according to the source code 
> that is defined as
>
> def square(a=1, name=''):
> """Make a square lattice."""
> return Monatomic(((a, 0), (0, a)), name=name)
>
> so I figured my definition should work for 3D.
>
> Could someone point out the obvious mistake I am making?
>
>
>
>

Re: [Kwant] forking smatrix.transmission

[Anton Akhmerov]

Hi Sergey,

What is the error traceback that you see? Is it anything related to
pickling? If so, then you should try the development version of Kwant,
which allows to pickle Kwant systems.

Best,
Anton

On Thu, Mar 30, 2017 at 10:41 PM, Sergey Slizovskiy  wrote:
> Dear all,
>
> I have tried a primitive forking of kwant computation with joblib package.
> This worked well for bands computation, here's the code, but it fails for
> smatrix.transmission.   I would be grateful for any hints and examples.
>
> Thank you,
>
> Sergey
>
> from joblib import Parallel, delayed
>
> def plot_bandstructure(flead, momenta, args, max=0.05):
> global  processors
> switch=False
> cutoff=0.001
> bands = kwant.physics.Bands(flead, args=args)
> energies = Parallel(n_jobs=processors)(delayed(bands)(k)  for k in
> momenta)
>
> pyplot.figure()
> pyplot.plot(momenta, energies)
> pyplot.ylim([-max,max])
> pyplot.xlabel("momentum [(lattice constant)^-1]")
> pyplot.ylabel("energy [t]")
> pyplot.show(block=False)
>
> def transm(sys, energy, lead1,lead2,args):
> print("Smatrix for energy ", energy)
> smatrix = delayed(kwant.smatrix)(sys, energy,args=args)
> return smatrix.transmission(lead1, lead2)
>
> def plot_conductance(sys, lead1,lead2, energies,args,label):
> # Compute transmission as a function of energy
>
> global processors
>
> data = Parallel(n_jobs=processors)(delayed(transm)(sys, energy,lead1,
> lead2, args=args)  for energy in energies)
> print(data)
> pyplot.figure()
> pyplot.plot(energies, data)
> pyplot.xlabel("energy [t]")
> pyplot.ylabel("conductance [e^2/h]")
> pyplot.title(label)
> return data
>

Re: [Kwant] How to reduce memory requirement

[Anton Akhmerov]

Dear Mingkai Li,

Firstly, 3D systems will take a lot of memory, and most likely you will not
be able to straightforwardly increase the system size by 50x without
switching to out-of-core or distributed computations. You could improve the
memory usage somewhat by also trying the metis or scotch orderings. In 2D
they greatly outperform PORD. However I don't think this will bring you
close to achieving what you plan.

The part that takes memory is storing the factors of the LU-decomposition.
For a nested dissection (the asymptotically optimal algorithm), the memory
requirement would be ~L^4 log L, with L the linear size of the system, and
it would require petabytes of storage. Even switching to the most
memory-efficient ordering and discarding unused LU-factors (so effectively
doing RGF), you still need to store a dense matrix with size equal to the
cross-section of your system, which already gives you a petabyte.

Best,
Anton

On Tue, Mar 14, 2017 at 2:14 AM, Li Mingkai  wrote:
> Dear all,
>
> I try to build a bulk transport system with periodic boundary condition by
> 40 * 40 * 40 sites. It requests over 2 GB memory during transmission
> calculation. I curious why it produces so large memory requirement. Since
> I'm going to expand the system up to 50 times, I want to know any way to
> reduce the memory requirement and increase calculation speed. I'm using
> kwant 1.1.2 and MUMPS with default PORD ordering. The system building code
> is
>
> syst = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((0, cube.NY,
> 0)), lat.vec((0, 0, cube.NZ
>
> for i in range(cube.NX):
> for j in range(cube.NY):
> for k in range(cube.NZ):
> syst[(lat(i, j, k))] = -cube.data[i, j, k] + 6 * t
>
> syst[lat.neighbors(1)] = -t
>
> lead = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((-1, 0, 0)),
> lat.vec((0, cube.NY, 0)), lat.vec((0, 0, cube.NZ
> lead[lat.shape(lead_shape, (0, 0, 0))] = 6 * t
> lead[lat.neighbors(1)] = -t
>
> syst=wraparound(syst)
> syst.attach_lead(wraparound(lead,keep=0))
> syst.attach_lead(wraparound(lead,keep=0).reversed())
>
> --
> Mingkai Li

Re: [Kwant] Problems with memory allocation (?) in KWANT 1.3

[Anton Akhmerov]

Hi Camilla,

It appears that for some reason you have an old version of tinyarray;
please also install the development version of tinyarray using pip3
install git+https://gitlab.kwant-project.org/kwant/tinyarray.git@master

Best,
Anton

On Fri, Mar 10, 2017 at 11:09 AM, Camilla Espedal
 wrote:
> Hi and thanks,
>
>
> I tried installing the newest Kwant version now, using:
>
> pip3 install git+https://gitlab.kwant-project.org/kwant/kwant.git@master
>
> but now I get a new error message (same code):
>
>
>   File "/home/camilla/Kwant/Tests/conservation.py", line 91, in 
> main()
>   File "/home/camilla/Kwant/Tests/conservation.py", line 68, in main
> syst = make_system(t=1.0, W=40, L=40, lead_W = 40)
>   File "/home/camilla/Kwant/Tests/conservation.py", line 56, in make_system
> syst.attach_lead(lead)
>   File "/home/camilla/.local/lib/python3.5/site-packages/kwant/builder.py",
> line 1402, in attach_lead
> self.leads.append(BuilderLead(lead_builder, tuple(interface)))
>   File "/home/camilla/.local/lib/python3.5/site-packages/kwant/builder.py",
> line 566, in __init__
> self.interface = tuple(sorted(interface))
> TypeError: unorderable types: tinyarray.ndarray_int() <
> tinyarray.ndarray_int()
> [Finished in 1.4s with exit code 1]
>
>
> should I have installed differently. According to pip3 I have installed
> kwant (1.3.0a0.dev123+g3961899)
>
>
> Sincerely,
>
> Camilla

Re: [Kwant] Kwant 1.3

[Anton Akhmerov]

Hi Ali,

Take a look at the new version of the new version of the
superconductor tutorial adapted for Kwant 1.3:
https://kwant-project.org/doc/dev/tutorial/tutorial5

Unfortunately, we don't have a tutorial for operators yet, it's in the
works. The best we have so far is just the documentation over here:
https://kwant-project.org/doc/dev/reference/kwant.operator

Best,
Anton

On Sat, Mar 4, 2017 at 12:35 PM, Ali Asgharpour
 wrote:
> Hello all,
>
> I was wondering if you'd mind explaining more about conservation law and
> operators you have added to Kwant 1.3. Would you please explain what
> conservation law does? Is it for just leads or it can be used for systems
> too? I have seen an example which stated in previous emails where you used
> conservation_law=-sz and time_reversal=s0. I want to know what happened if I
> use sx instead of sz numerically.
>
> In addition, I want to know more about Current and Source operators. It
> would be great if you implement them in a small code like conservation law
> one as an example.
>
> Your time and consideration are greatly appreciated.
>
> Best regards,
>
> Ali

Re: [Kwant] Smatrix in hall bar

[Anton Akhmerov]

Hi Ali,

We subtract the sums because that is the definition of the conductance
matrix. However this complete function really should be written as
kwant.smatrix(syst, args=[p]).conductance_matrix() (see also the
documentation of the conductance_matrix from Kwant).

Best,
Anton

On Wed, Feb 15, 2017 at 3:43 PM, Ali Asgharpour
 wrote:
> Hello all,
>
> I have a question about one of the functions in this link:
> https://nbviewer.jupyter.org/url/topocondmat.org/notebooks/w3_pump_QHE/Laughlinargument.ipynb
>
> For plotting Hall conductance and longitudinal one, there is a section where
> you defined a G function, which I put it below:
>
> def G(syst, p):
> smatrix = kwant.smatrix(syst, args=[p])
> G = [[smatrix.transmission(i, j) for i in range(num_leads)] for j in
> range(num_leads)]
> G -= np.diag(np.sum(G, axis=0))
> return calculate_sigmas(G)
>
>
> I want to know why we need this line " G -= np.diag(np.sum(G, axis=0)) "
> since G matrix is already calculated in the previous line?
>
> Your time and consideration are greatly appreciated.
>
> Best regards,
>
> Ali
>
>
>
>
>
>
>

Re: [Kwant] Regarding smatrix and spin

[Anton Akhmerov]

Hi Camilla,

The reason you are seeing different scattering matrices is because
without knowledge of the conservation laws Kwant orders the modes by
their velocities first, and not by their internal degrees of freedom.
Transmissions should be the same, however.

The updated way to specify the conservation laws is much simpler than
in the first draft implementation that we shared. You just need to
provide the conserved quantity on each site when constructing the
system (e.g. s_z). This operator must have integer eigenvalues. Then
the blocks correspond to the eigenvalues of that operator sorted in
ascending order (so (0, 1) comes first, (1, 0) second).

Best,
Anton


On Tue, Feb 14, 2017 at 11:26 AM, Camilla Espedal
<camilla.espe...@ntnu.no> wrote:
> Hi Anton,
>
> Thank you! It is working now. I have one question though, about the 
> conservation laws and Kwant.
>
> So, I tried to add a really small zeeman-field in the leads (delta_z * 
> sigma_z) sot that the two spin states are not degenerate. The idea was that 
> it would be small enough to not affect the physics, but still force the 
> spin-basis. However when I compare the smatrix.submatrix(1,0) in the case 
> where I add the small Zeeman field, and the one I get using the conservation 
> laws, they are not the same.
>
> Do you have any idea as to why?
>
> Also, what is the convention in the numbering of the spins (q_0 and q_1)? Is 
> q_0 = q_1 = 0 upup? It seems that the syntax has changed since you used 
> set_symmetry before?
>
> Best,
> Camilla
> 
> Fra: Anton Akhmerov <anton.akhme...@gmail.com>
> Sendt: 14. februar 2017 09:48
> Til: Camilla Espedal
> Kopi: kwant-discuss@kwant-project.org
> Emne: Re: [Kwant] Regarding smatrix and spin
>
> Hi Camilla,
>
> Sure. You need to install a development version of tinyarray. If you
> got everything via conda, then you can do " conda install -c kwant
> 'tinyarray==dev' ". If you installed kwant via pip, do " pip install
> git+https://gitlab.kwant-project.org/kwant/tinyarray.git@master "
>
> Let me know if it works (both the installation and using the conservation 
> laws).
>
> Best,
> Anton
>
> On Tue, Feb 14, 2017 at 9:40 AM, Camilla Espedal
> <camilla.espe...@ntnu.no> wrote:
>> Dear Anton,
>>
>> sorry for troubling you again.
>>
>> so I finally got hold of a linux-computer, and managed to run the commands
>> and install all the things in the notebook. When I try to run the code you
>> have in your notebook however, I get the following error message:
>>
>>   File "/home/camilla/.local/lib/python3.5/site-packages/kwant/builder.py",
>> line 1371, in attach_lead
>> self.leads.append(BuilderLead(lead_builder, tuple(interface)))
>>   File "/home/camilla/.local/lib/python3.5/site-packages/kwant/builder.py",
>> line 565, in __init__
>> self.interface = tuple(sorted(interface))
>> TypeError: unorderable types: tinyarray.ndarray_int() <
>> tinyarray.ndarray_int()
>>
>> Do you know what the problem could be?
>>
>> Best,
>>
>> Camilla
>>
>> The code if needed:
>>
>> import kwant
>> import tinyarray as ta
>> import numpy as np
>> from scipy import sparse
>> from matplotlib import pyplot
>> import matplotlib
>>
>> s0 = ta.array([[1, 0], [0, 1]])
>> sx = ta.array([[0, 1], [1, 0]])
>> sy = ta.array([[0, 1j], [-1j, 0]])
>> sz = ta.array([[1, 0], [0, -1]])
>>
>>
>>
>> # Adapted from https://kwant-project.org/doc/1/tutorial/tutorial2
>>
>> def make_system(t=1.0, W=10, L=10):
>> # Now we must specify the number of orbitals per site.
>> lat = kwant.lattice.square(norbs=2)
>> syst = kwant.Builder()
>>
>> syst[(lat(x, y) for x in range(L) for y in range(W))] = \
>> lambda s, alpha, E_z: 4 * t * s0 + E_z * sz
>> syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
>> lambda s1, s2, alpha, E_z: -t * s0 - 1j * alpha * sy
>> syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
>> lambda s1, s2, alpha, E_z: -t * s0 + 1j * alpha * sx
>>
>> # The new bit: specifying the conservation law.
>> lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)),
>>  conservation_law=-sz, time_reversal=s0)
>> lead[(lat(0, j) for j in range(W))] = 4 * t * s0
>> lead[lat.neighbors()] = -t * s0 # Note: no spin-orbit in the lead.
>>
>> syst.attach_lead(lead)
>> syst.attach_lead(lead.reversed())
>>
>> syst = syst.finalized()
>>
&g

Re: [Kwant] Regarding smatrix and spin

[Anton Akhmerov]

Hi Camilla,

Sure. You need to install a development version of tinyarray. If you
got everything via conda, then you can do " conda install -c kwant
'tinyarray==dev' ". If you installed kwant via pip, do " pip install
git+https://gitlab.kwant-project.org/kwant/tinyarray.git@master "

Let me know if it works (both the installation and using the conservation laws).

Best,
Anton

On Tue, Feb 14, 2017 at 9:40 AM, Camilla Espedal
<camilla.espe...@ntnu.no> wrote:
> Dear Anton,
>
> sorry for troubling you again.
>
> so I finally got hold of a linux-computer, and managed to run the commands
> and install all the things in the notebook. When I try to run the code you
> have in your notebook however, I get the following error message:
>
>   File "/home/camilla/.local/lib/python3.5/site-packages/kwant/builder.py",
> line 1371, in attach_lead
> self.leads.append(BuilderLead(lead_builder, tuple(interface)))
>   File "/home/camilla/.local/lib/python3.5/site-packages/kwant/builder.py",
> line 565, in __init__
> self.interface = tuple(sorted(interface))
> TypeError: unorderable types: tinyarray.ndarray_int() <
> tinyarray.ndarray_int()
>
> Do you know what the problem could be?
>
> Best,
>
> Camilla
>
> The code if needed:
>
> import kwant
> import tinyarray as ta
> import numpy as np
> from scipy import sparse
> from matplotlib import pyplot
> import matplotlib
>
> s0 = ta.array([[1, 0], [0, 1]])
> sx = ta.array([[0, 1], [1, 0]])
> sy = ta.array([[0, 1j], [-1j, 0]])
> sz = ta.array([[1, 0], [0, -1]])
>
>
>
> # Adapted from https://kwant-project.org/doc/1/tutorial/tutorial2
>
> def make_system(t=1.0, W=10, L=10):
> # Now we must specify the number of orbitals per site.
> lat = kwant.lattice.square(norbs=2)
> syst = kwant.Builder()
>
> syst[(lat(x, y) for x in range(L) for y in range(W))] = \
> lambda s, alpha, E_z: 4 * t * s0 + E_z * sz
> syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
> lambda s1, s2, alpha, E_z: -t * s0 - 1j * alpha * sy
> syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
> lambda s1, s2, alpha, E_z: -t * s0 + 1j * alpha * sx
>
> # The new bit: specifying the conservation law.
> lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)),
>  conservation_law=-sz, time_reversal=s0)
> lead[(lat(0, j) for j in range(W))] = 4 * t * s0
> lead[lat.neighbors()] = -t * s0 # Note: no spin-orbit in the lead.
>
> syst.attach_lead(lead)
> syst.attach_lead(lead.reversed())
>
> syst = syst.finalized()
>
>
> return syst
>
> syst = make_system(t=1.0, W=10, L=10)
> energies = np.linspace(0, 1, 200)
> smatrices = [kwant.smatrix(syst, energy, args=(0.2, 0.05)) for energy in
> energies]
>
> fig = pyplot.figure(figsize=(13, 8))
> ax = fig.add_subplot(1, 1, 1)
>
> # Like previously smatrix.transmission(lead1, lead0) is transmission from
> lead0 to lead1
> ax.plot(energies, [smatrix.transmission(1, 0) for smatrix in smatrices],
> label='total')
>
> # The new bit: smatrix.transmission((lead1, q1), (lead0, q0)) is the
> transmission from the
> # q0 block of the lead0 into the q1 block of lead1. The subblock ordering is
> same as we used
> # in set_symmetry.
> ax.plot(energies, [smatrix.transmission((1, 0), (0, 0)) for smatrix in
> smatrices], label='$G_{↑↑}$')
> ax.plot(energies, [smatrix.transmission((1, 1), (0, 0)) for smatrix in
> smatrices], label='$G_{↑↓}$')
> ax.plot(energies, [smatrix.transmission((1, 0), (0, 1)) for smatrix in
> smatrices], label='$G_{↓↑}$')
> ax.plot(energies, [smatrix.transmission((1, 1), (0, 1)) for smatrix in
> smatrices], label='$G_{↓↓}$')
> ax.set_ylabel('$G [e^2/h]$', fontsize='xx-large')
> ax.set_xlabel('$E/t$', fontsize='xx-large')
> ax.legend(fontsize='x-large');
>
>
> On 24. jan. 2017 13:18, Anton Akhmerov wrote:
>
> OK, please double-check the remaining simulation parameters.
>
> Anton
>
> On Tue, Jan 24, 2017 at 1:00 PM, Camilla Espedal
> <camilla.espe...@ntnu.no> wrote:
>
> Dear Anton,
>
>
>
> I triend changing it, but that does not solve the problem.
>
>
>
> Best,
>
> Camilla
>
> From: Anton Akhmerov [mailto:anton.akhmerov...@gmail.com]
> Sent: 24. januar 2017 12:36
>
>
> To: Camilla Espedal <camilla.espe...@ntnu.no>
> Cc: kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] Regarding smatrix and spin
>
>
>
> Dear Camilla,
>
>
>
> Could the difference originate from you using a lattice constant of 2
> instead of 1?
>
>
>
> Anton
>
>
>
> On Tue, Jan 24, 2017, 10:41 C

Re: [Kwant] Regarding tranlational invariance along 2 directions

[Anton Akhmerov]

Dear Anant,

You are searching for the wraparound module built by Christoph and
previously announce on this mailing list. Please see the original
announcement over here:
https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00745.html

Once again I would like to remind you to try searching for questions
on this mailing list before asking questions.

Anton

On Tue, Jan 24, 2017 at 8:04 AM, ANANT VIJAY  wrote:
> Hello everyone,
>I am writing this to ask about how can I finalize the
> sytem with translational invariance along two directions.
> Code that I am using is the following as I want to make a slab system for
> band structure.
>
> sym = kwant.TranslationalSymmetry([1, 0, 0], [0, 1, 0])
>
> Error:
>
> ValueError: Currently, only builders without or
> with a 1D translational
>  symmetry can be finalized.
>
> Is there any solution to this?
>
> I also happen to get one code from notes on internet by
> edx course where thwe following code has been used using wraparound
> function.
>
> I did not get much out of it. Also, I checked I don't have this function in
> my kwant module. how can I get this functions module?
> Please anyone can explian what this function is doing in the code:
>
> def make_scatter_sys():
> from functions import wraparound
> syst = wraparound(bhz(Z=1, system_type='slab'))
> syst.attach_lead(make_lead())
> syst.attach_lead(wraparound(bhz(), keep=0))
> syst = syst.finalized()
> syst.leads[0] = TRIInfiniteSystem(syst.leads[0], trs)
> return syst
>
>
> ANANT VIJAY VARMA
> Research Scholar
> IISER-KOLKATA
> WEST BENGAL.

Re: [Kwant] Regarding smatrix and spin

[Anton Akhmerov]

Dear Camilla,

Could the difference originate from you using a lattice constant of 2
instead of 1?

Anton

On Tue, Jan 24, 2017, 10:41 Camilla Espedal <camilla.espe...@ntnu.no> wrote:

> Dear Anton,
>
> Thanks again for all your help. I will try to do it the linux way. Just
> one more thing regarding this. I wrote a script in the old Kwant to find
> the total conductance and compare it to the one in your notebook. While the
> two plots are qualitatively similar, they are not the same. Am I missing
> something, or am I calculating different things?
>
> Best, Camilla
>
> (my code):
>
> # Tutorial 2.3.1. Matrix structure of on-site and hopping elements
> # 
> #
> # Physics background
> # --
> #  Gaps in quantum wires with spin-orbit coupling and Zeeman splititng,
> #  as theoretically predicted in
> #   http://prl.aps.org/abstract/PRL/v90/i25/e256601
> #  and (supposedly) experimentally oberved in
> #   http://www.nature.com/nphys/journal/v6/n5/abs/nphys1626.html
> #
> # Kwant features highlighted
> # --
> #  - Numpy matrices as values in Builder
>
> import kwant
>
> # For plotting
> import matplotlib.pyplot as plt
>
> # For matrix support
> import tinyarray
> import numpy as np
>
> # define Pauli-matrices for convenience
> sigma_0 = tinyarray.array([[1, 0], [0, 1]])
> sigma_x = tinyarray.array([[0, 1], [1, 0]])
> sigma_y = tinyarray.array([[0, 1j], [-1j, 0]])
> sigma_z = tinyarray.array([[1, 0], [0, -1]])
>
>
> def make_system(a=2, t=1.0, alpha=0.1, e_z=0.05, W=10, L=10):
> # Start with an empty tight-binding system and a single square lattice.
> # `a` is the lattice constant (by default set to 1 for simplicity).
> lat = kwant.lattice.square(a)
>
> sys = kwant.Builder()
>
>  Define the scattering region. 
> sys[(lat(x, y) for x in range(L) for y in range(W))] = \
> 4 * t * sigma_0 + e_z * sigma_z
> # hoppings in x-direction
> sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
> -t * sigma_0 - 1j * alpha * sigma_y
> # hoppings in y-directions
> sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
> -t * sigma_0 + 1j * alpha * sigma_x
>
>  Define the left lead. 
> lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
>
> lead[(lat(0, j) for j in range(W))] = 4 * t * sigma_0
> # hoppings in x-direction
> lead[lat.neighbors()] = \
> -t * sigma_0
>
>  Attach the leads and return the finalized system. 
> sys.attach_lead(lead)
> sys.attach_lead(lead.reversed())
>
> return sys
>
>
> def plot_conductance(sys, energies):
> # Compute conductance
> data = []
> for energy in energies:
> smatrix = kwant.smatrix(sys, energy)
> data.append(smatrix.transmission(1, 0))
>
> pyplot.figure()
> pyplot.plot(energies, data)
> pyplot.xlabel("energy [t]")
> pyplot.ylabel("conductance [e^2/h]")
> pyplot.show()
>
>
> def main():
> sys = make_system()
>
> # Check that the system looks as intended.
> kwant.plot(sys)
>
> # Finalize the system.
> sys = sys.finalized()
> energies = np.linspace(0, 1, 200)
> smatrices = [kwant.smatrix(sys, energy) for energy in energies]
>
> fig = plt.figure(figsize=(13, 8))
> ax = fig.add_subplot(1, 1, 1)
>
> ax.plot(energies, [smatrix.transmission(1,0) for smatrix in
> smatrices], label='total')
>
> ax.set_ylabel('$G [e^2/h]$', fontsize='xx-large')
> ax.set_xlabel('$E/t$', fontsize='xx-large')
> ax.legend(fontsize='x-large')
>
> plt.show()
>
>
>
> # Call the main function if the script gets executed (as opposed to
> imported).
> # See <http://docs.python.org/library/__main__.html>.
> if __name__ == '__main__':
> main()
>
> -Original Message-
> From: anton.akhme...@gmail.com [mailto:anton.akhme...@gmail.com] On
> Behalf Of Anton Akhmerov
> Sent: 17. januar 2017 10:48
> To: Camilla Espedal <camilla.espe...@ntnu.no>
> Cc: kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] Regarding smatrix and spin
>
> Dear Camilla,
>
> It seems that you are trying to install Kwant on windows. This is a very
> hard task, and I fear none of the Kwant developers has enough knowledge of
> it right now (our Windows packages are built by Christoph Gohlke, see [1]
> for the build environment description). However if you are using windows
> 10, I suggest to try to install Kwant using the windows subsystem for
> linux. That way th

Re: [Kwant] Units of the wavefunction and velocity

[Anton Akhmerov]

Hi Harshad,

> I was able to obtain the velocities from the modes object. I confirmed
> that modes.velocities is just the inverse of integral |휓|2 over the unit
> cell. If I did the math correctly, I got the units of modes.velocities to be
> eV (energy units of my Hamiltonian). How can I convert this to lets say m/s?

Sorry that I forgot to clarify that the unit of length in this case is
the side of the lead unit cell. Multiplying by that size should allow
to easily translate your velocity into m/s.

Anton

[Kwant] Installation instructions for Kwant on WSL

[Anton Akhmerov]

Hi everyone,

Andrey Antipov kindly shared his installation instructions for Kwant
using the Windows subsystem for Linux. If you are interested, please
check it out over here:
https://gist.github.com/aeantipov/3839d09529cb5dd2cedc00b76ccd2eb8

Best,
Anton

Re: [Kwant] Units of the wavefunction and velocity

[Anton Akhmerov]

> Sorry I made a mistake in the units of 'I'. Isn't 'I' here the probability 
> current and not the charge current? We would just have to multiply 'I' by e 
> to get the charge current, right?

Yes, that is correct. The equation that you quoted is indeed the
probability current.

Re: [Kwant] Units of the wavefunction and velocity

[Anton Akhmerov]

Dear Hardshad,

I might be missing some factors here, but the unit of current is elementary
charge / unit of time. The unit of time is hbar / unit of energy that you
used in defining your tight-binding model.

Best,
Anton

On Mon, Jan 23, 2017 at 3:49 PM, Harshad Sahasrabudhe 
wrote:

> Hi,
>
> I was wondering what units would the wavefunction obtained from Kwant
> have?  I was thinking they would have the units 1/sqrt(nm.eV) (since my
> energies are in eV and lengths are in nm) as the modes are normalized
> according to
>
> [image: Inline image 1]
>
> How would I calculate the velocity of the mode, if the modes are
> normalized to carry unit current?
>
> Thanks,
> Harshad
>

Re: [Kwant] Regarding smatrix and spin

[Anton Akhmerov]

Dear Camilla,

It seems that you are trying to install Kwant on windows. This is a
very hard task, and I fear none of the Kwant developers has enough
knowledge of it right now (our Windows packages are built by Christoph
Gohlke, see [1] for the build environment description). However if you
are using windows 10, I suggest to try to install Kwant using the
windows subsystem for linux. That way the standard Ubuntu build
procedure should work for you.

Best,
Anton

[1]: http://www.lfd.uci.edu/~gohlke/pythonlibs/

On Mon, Jan 16, 2017 at 9:45 AM, Camilla Espedal
<camilla.espe...@ntnu.no> wrote:
> Thanks a lot. I tried to install the cons_laws_combined, but I get the 
> following error message:
>
> "LINK: fatal error LNK1181: cannot open input file 'lapack.lib'"
>
> Is there some package or installation I am missing?
>
> Best regards,
> Camilla
>
> -Original Message-
> From: anton.akhme...@gmail.com [mailto:anton.akhme...@gmail.com] On Behalf Of 
> Anton Akhmerov
> Sent: 8. januar 2017 16:35
> To: Tómas Örn Rosdahl <torosd...@gmail.com>
> Cc: Camilla Espedal <camilla.espe...@ntnu.no>; kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] Regarding smatrix and spin
>
> Hi Camilla, everyone,
>
> I've slightly modified Tómas's example to a case where the spins do get 
> coupled, check it out:
> http://nbviewer.jupyter.org/url/antonakhmerov.org/misc/spin_conductance.ipynb
>
> I've also provided more detailed installation instructions in the notebook.
>
> Cheers,
> Anton
>
> On Sun, Jan 8, 2017 at 2:45 PM, Tómas Örn Rosdahl <torosd...@gmail.com> wrote:
>> Dear Camilla,
>>
>> For a Hamiltonian with degeneracies due to a conservation law, the
>> scattering states will in general not have a definite value of the
>> conservation law. In your case, Kwant returns scattering states that
>> are arbitrary linear combinations of spin up and down, so it is not
>> possible to label the amplitudes in the scattering matrix by spin.
>>
>> However, in Kwant 1.3 a feature will be added that allows for the
>> construction of scattering states with definite values of a
>> conservation law. See here for an explanation of the basic idea behind the 
>> algorithm.
>>
>> We're currently working on implementing this feature in Kwant itself.
>> The good news is that we're practically done - here is a link to a git
>> repo with a functioning implementation. After you clone the repo,
>> check out the branch cons_laws_combined, which contains a version of
>> Kwant with conservation laws implemented. This notebook contains a
>> simple example to illustrate how to work with conservation laws and the 
>> scattering matrix.
>>
>> I invite you and anyone else who is interested to give it a try. We'd
>> appreciate any feedback!
>>
>> In your case specifically, there would be two projectors in the new
>> implementation - P0 which projects out the spin up block, and P1 that
>> projects out the spin down block. If they are specified in this order,
>> then the spin up and down blocks in the Hamiltonian have block indices
>> 0 and 1, respectively. In the new implementation, it is possible to
>> ask for subblocks of the scattering matrix relating not only any two
>> leads, but also any two conservation law blocks in any leads. To get
>> the reflection amplitude of an incident spin up electron from lead 0
>> into an outgoing spin down electron in lead 0, you could simply do
>> smat.submatrix((0, 1), (0, 0)). Here, the arguments are tuples of indices 
>> (lead index, block index).
>>
>> Best regards,
>> Tómas
>>
>> On Fri, Jan 6, 2017 at 3:46 PM, Camilla Espedal
>> <camilla.espe...@ntnu.no>
>> wrote:
>>>
>>> Hi again,
>>>
>>>
>>>
>>> This question is basically the same as this:
>>> https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00076
>>> .html
>>>
>>>
>>>
>>> I want to calculate some things using the scattering matrix. I
>>> started out with a very simple system, most basic two-terminal
>>> system. For some energy there is one propagating mode. I now add
>>> matrix structure to the mix (just multiply by s_0 everywhere) and
>>> there are now 2 propagating modes (which makes sense).
>>>
>>>
>>>
>>> Now, if I look at the reflection coefficients for lead 0 by using
>>> submatrix(0,0), it is now a 2x2 matrix after I introduced the
>>> matrices. How are the elements ordered? Is it
>>>
>>>
>>>
>>> [[r_upup, r_updown],[r_downup, r_downdown]]
>>>
>>>
>>>
>>> I know that I could make two lattices, but since I do not plan to use
>>> the other functions such as transmission. I  just want the smatrix.
>>>
>>>
>>>
>>> Hope you can help me, and thanks in advance.
>>>
>>>
>>>
>>> Best regards,
>>>
>>> Camilla
>>
>>

Re: [Kwant] the path in the Brillouin zone

[Anton Akhmerov]

Dear Weiyuan,

Since there's only a single translationally invariant direction, the
Brillouin zone is one-dimensional, and therefore you are viewing the
complete 1D Brillouin zone. If there was a finite order translation
invariance in the transverse direction like in a nanotube, each band
would correspond to a different discrete eigenvalue of momentum, as
described in a lot of educational materials, e.g. [1]. However in a
nanoribbon the momentum in a transverse direction is not conserved,
and therefore the 1D bands cannot in general be mapped on a 1D path in
a 2D Brillouin zone.

Best,
Anton

[1]: https://nanoelectronics.unibas.ch/education/Nanotubes/LCAO-NT.pdf

On Tue, Jan 17, 2017 at 5:59 AM, Weiyuan Tong
 wrote:
> Dear all,
> I used kwant to calculate the band structure of zigzag and armchair graphene
> nanoribbons. For zigzag graphene nanoribbons, we have sym=
> kwant.TranslationalSymmetry(lat.a.vec((-1, 0))). For armchair graphene
> nanoribbons, we have sym= kwant.TranslationalSymmetry(lat.a.vec((-1, 2))).
> My question is: which is the path in the Brillouin zone that we use and
> label the x-axis accordingly? My figure is attached. Thanks in advance!
> Weiyuan

Re: [Kwant] Regarding smatrix and spin

[Anton Akhmerov]

Hi Camilla, everyone,

I've slightly modified Tómas's example to a case where the spins do
get coupled, check it out:
http://nbviewer.jupyter.org/url/antonakhmerov.org/misc/spin_conductance.ipynb

I've also provided more detailed installation instructions in the notebook.

Cheers,
Anton

On Sun, Jan 8, 2017 at 2:45 PM, Tómas Örn Rosdahl  wrote:
> Dear Camilla,
>
> For a Hamiltonian with degeneracies due to a conservation law, the
> scattering states will in general not have a definite value of the
> conservation law. In your case, Kwant returns scattering states that are
> arbitrary linear combinations of spin up and down, so it is not possible to
> label the amplitudes in the scattering matrix by spin.
>
> However, in Kwant 1.3 a feature will be added that allows for the
> construction of scattering states with definite values of a conservation
> law. See here for an explanation of the basic idea behind the algorithm.
>
> We're currently working on implementing this feature in Kwant itself. The
> good news is that we're practically done - here is a link to a git repo with
> a functioning implementation. After you clone the repo, check out the branch
> cons_laws_combined, which contains a version of Kwant with conservation laws
> implemented. This notebook contains a simple example to illustrate how to
> work with conservation laws and the scattering matrix.
>
> I invite you and anyone else who is interested to give it a try. We'd
> appreciate any feedback!
>
> In your case specifically, there would be two projectors in the new
> implementation - P0 which projects out the spin up block, and P1 that
> projects out the spin down block. If they are specified in this order, then
> the spin up and down blocks in the Hamiltonian have block indices 0 and 1,
> respectively. In the new implementation, it is possible to ask for subblocks
> of the scattering matrix relating not only any two leads, but also any two
> conservation law blocks in any leads. To get the reflection amplitude of an
> incident spin up electron from lead 0 into an outgoing spin down electron in
> lead 0, you could simply do smat.submatrix((0, 1), (0, 0)). Here, the
> arguments are tuples of indices (lead index, block index).
>
> Best regards,
> Tómas
>
> On Fri, Jan 6, 2017 at 3:46 PM, Camilla Espedal 
> wrote:
>>
>> Hi again,
>>
>>
>>
>> This question is basically the same as this:
>> https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00076.html
>>
>>
>>
>> I want to calculate some things using the scattering matrix. I started out
>> with a very simple system, most basic two-terminal system. For some energy
>> there is one propagating mode. I now add matrix structure to the mix (just
>> multiply by s_0 everywhere) and there are now 2 propagating modes (which
>> makes sense).
>>
>>
>>
>> Now, if I look at the reflection coefficients for lead 0 by using
>> submatrix(0,0), it is now a 2x2 matrix after I introduced the matrices. How
>> are the elements ordered? Is it
>>
>>
>>
>> [[r_upup, r_updown],[r_downup, r_downdown]]
>>
>>
>>
>> I know that I could make two lattices, but since I do not plan to use the
>> other functions such as transmission. I  just want the smatrix.
>>
>>
>>
>> Hope you can help me, and thanks in advance.
>>
>>
>>
>> Best regards,
>>
>> Camilla
>
>

Re: [Kwant] STM tip as lead

[Anton Akhmerov]

Hi Eric,

Do you mind providing more information about what you are trying to
achieve? An STM tip is much simpler than an actual lead, so one could
avoid the complexity potentially.

Best,
Anton

On Sun, Jan 8, 2017 at 1:28 AM, Eric Mascot  wrote:
> Hello
>
> I was wondering, is it possible to have a lead that hovers above a 2D
> system?
> I'm guessing that I have to use a 3D system where the system is in the XY
> plane and the lead has a translational symmetry in the Z direction.
>
> Thanks,
> Eric

Re: [Kwant] problem with installing kwant on supercomputer

[Anton Akhmerov]

Dear Shizeng,

The error messages indicate that the compilation didn't work
successfully. While I'm not entirely sure what goes wrong, I'd like to
suggest a simpler alternative. If you are using anaconda on a *nix
system, the most straightforward way to obtain kwant is to use kwant
provided via the community-maintained conda-forge channel. Just
execute

conda install -c conda-forge kwant

and you are done.

Best,
Anton

On Tue, Jan 3, 2017 at 10:51 PM, Shizeng Lin  wrote:
> Hi,
> I try to install kwant on a supercomputer in my institute. I do not
> have root access, so I built the package from source.
> I am using python/3.5-anaconda-4.1.1.
> I installed tinyarray and kwant as a user module, with the following command
>
> pip install --user tinyarray/
> pip install --user kwant/
>
> The installation completed successfully.
>
> But when I try to run the test
>
> python -c 'import kwant; kwant.test()'
>
> I got the following error message:
>
> Error message=
> Traceback (most recent call last):
>
>   File "/turquoise/users/shizeng/kwant/kwant/__init__.py", line 14, in 
> 
>
> from . import _system
>
> ImportError: cannot import name '_system'
>
>
> During handling of the above exception, another exception occurred:
>
>
> Traceback (most recent call last):
>
>   File "", line 1, in 
>
>   File "/turquoise/users/shizeng/kwant/kwant/__init__.py", line 22, in 
> 
>
> raise ImportError(msg)
>
> ImportError: Error importing Kwant:
>
> You should not try to import Kwant from its source directory.
>
> Please exit the Kwant source distribution directory, and relaunch
>
> your Python intepreter from there.
>
> ==End=
> I followed the instruction to launch the test at different directory.
> but still got the same error message. I am appreciated if someone can
> help me to figure out the problem.
> Thanks.
> Shizeng

Re: [Kwant] On-site potential only at the edge of the sample in Kwant

[Anton Akhmerov]

Dear Dongwook,

It is the easiest to achieve what you want while you are defining your
scattering region, specifically by using the Builder.degree method. So
if you have a builder syst with all the hoppings already added, your
desired result would be achieved by
syst[(site for site in syst.sites() if syst.degree(site))) != 4)] =
boundary_potential

Best,
Anton

On Tue, Jan 3, 2017 at 4:08 AM, Dongwook Go  wrote:
> Dear developers,
>
> I’m interested in investigating the influence of the edge potential of the 
> sample.
> For example, in an example by tutorial “quantum_wire.py”, I would like to put 
> on-site potential only at the edge atoms.
>
> One could determine whether the site is edge or not by reading the number of 
> neighbors of the site by using “degree”.
> For example, if the degree is 3 there is an on-site potentials , else 0.
> Is there any example code from which I can learn to use this function?
>
>
> Best,
> Dongwook.

Re: [Kwant] Transmission in ZGNR

[Anton Akhmerov]

Dear Abhishek Sharan,

Please see the previous answers to similar questions:
https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00983.html

Best regards,
Anton Akhmerov

On Sat, Dec 3, 2016 at 2:35 AM, Abhishek Sharan <abhksha...@gmail.com> wrote:
> Hi!
>
> I am new to using kwant, so some of my questions might be naive.
>
> I need to compute conductance (or transmission) and fano factor of the shot
> noise of zig-zag graphene nanoribbon for a particular geometry and a set
> gate voltage and a set fermi energy.
>
> I can calculate the conductance of ZGNR for a particular geometry as a
> function of energy.
> I need to understand two things:
>
> 1. I assume setting the gate voltage means fixing the energy of the
> electrons and finding the transmission at that energy (right?)
>
> 2. I do not understand how do I set the fermi level in the system?
>
> I would appreciate if someone could clarify on this.
>
> --
> Regards
>

Re: [Kwant] Question regarding orbitals and scattering matrix

[Anton Akhmerov]

Hi Camilla,

The scattering matrix contains the transmission/reflection amplitudes
of propagating modes. Those are in general superpositions of all the
orbitals present in the system. This is why the probability to
scattering from orbital 1 in one lead to orbital 2 in another lead
isn't defined, since it will change depending on e.g. distance from
the scattering region.

Hope that clarifies things,
Anton

On Fri, Nov 25, 2016 at 2:48 AM, Camilla Espedal
 wrote:
> Hi again,
>
>
>
> This question is sort of related to one of mye previous questions
> (https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01029.html)
> but I am writing it as a new one. The questions relates to the sites in the
> lead unit cell and the scattering matrix.
>
>
>
> Say I have a system (simple 2D wire), where I have 3 sites in my unit cell
> (3 orbitals). I then make a system at an energy which is so low that I only
> have one propagating mode in the system. If I then run
>
>
>
> smatrix.lead_info[0]
>
>
>
> it returns a 2x3 matrix (which corresponds to left/right going and number of
> orbitals). If I then run
>
>
>
> smatrix.submatrix(0,0)
>
>
>
> I only get one element, which orbital does this element “belong” to, or how
> does it relate to these orbitals? Is the scattering matrix the same for all
> orbitals? Is the probability of scattering from say orbital 1 in one lead to
> orbital 2 in another the same as from orbital 1 to orbital 1?
>
>
>
> Hope my question makes sens.
>
>
>
> Best regards,
>
> Camilla

Re: [Kwant] smatrix for up and down spin channels in graphene with KM Hamiltonian

[Anton Akhmerov]

Hi Sudin,

Since transmission is quantized in one case and isn't quantized in the
other, it seems that in case II the Hamiltonian of the leads is
different from the Hamiltonian of the scattering region. Moreover I
see that in case I you have a term that is also present in the leads
and that is proportional to sigma_x. This term cannot be added to case
II since there it would be a hopping term between different leads (and
that's not allowed by Kwant). Physically, in case I you don't have a
conservation law in the leads, so you cannot separate the lead
Hamiltonian into that of different spins.

Hope that helps,
Anton

On Wed, Nov 23, 2016 at 6:06 AM, Sudin Ganguly  wrote:
> Dear Sir,
>
> Now I found the solution, I should have named the different kinds of
> lattices, like
>
> 
> lat_u = kwant.lattice.honeycomb(1, name='up')
> a,b = lat_u.sublattices
>
> lat_d = kwant.lattice.honeycomb(1, name='down')
> c,d = lat_d.sublattices
> ===
>
>
> Now I have another problem.
>
> I was trying to build up the Kane-Mele Hamiltonian on graphene.
>
> Case I: Followed by the  example
> "https://kwant-project.org/doc/1.0/tutorial/tutorial2;, I've included
> Rashba term, and got the two terminal conductance.
>
> Case II: Now I want the transmission coefficients from different spin
> species. So I build up the Hamiltonian for up and down spins.
>
> If I want to compare case I and case II, one needs to calculate
>
> for case I   smatrix(1,0)
>  and
> for case II smatrix(2,0) + smatrix(3,0) + smatrix(2,1) + smatrix(3,1).
>
> But these two are not coming same, which should be. When Rashba strength
> is zero, It's ok.
>
> I am not getting any clue why these two results are not coming out the same.
>
> I've appended two codes for case I and case II.
>
>
>
> With regards,
> Sudin
>
>
>
>
>
>
>
>> Dear sir,
>>
>> I've implemented the same on square lattice, and the conductance is coming
>> out as expected. But when I did the same on graphene, the expected result
>> is not coming.
>>
>> With Regards,
>> Sudin
>>
>>
>
>
> --
> 
> Sudin Ganguly
> Research Scholar
> Dept. of Physics
> IIT Guwahati
> Assam,India-781039
> ===

Re: [Kwant] Hopping between different lattices and families

[Anton Akhmerov]

Hi Antonio,

It's not really about Kwant, but rather how you define your system.
Specifically you have forgotten one pairing term. Using

sys[((lat_e(x,y), lat_h(x+1,y)) for x in range(L-1) for y in range(W))] =
Delta
sys[((lat_h(x,y), lat_e(x+1,y)) for x in range(L-1) for y in range(W))] =
-Delta

fixes the problem.

Best,
Anton

On Thu, Nov 3, 2016 at 9:52 AM, <antonio...@usp.br> wrote:

> Hi Anton,
>
> I want to implement a p-wave superconductor using two lattices (for
> electrons and holes). I've done that using Pauli matrices with no problem.
> The thing is that I expect to get Majorana zero energy modes, but  I don't.
> I belive that the problem is the way that I am defining the pairing, which
> does not take account for the lat_h(0,0) and lat_e(L,0) positions. I wonder
> if any of you have some tips to implement this coupling with this approach
> (two lattices) to get the zero energy modes.
>
> Best,
> Antônio.
>
> --
>
> *De: *"Anton Akhmerov" <anton.akhme...@gmail.com>
> *Para: *antonio...@usp.br
> *Cc: *kwant-discuss@kwant-project.org, "Bas Nijholt" <basnijh...@gmail.com>,
> "Camilla Espedal" <camilla.espe...@ntnu.no>
> *Enviadas: *Quarta-feira, 2 de Novembro de 2016 22:26:30
>
> *Assunto: *Re: [Kwant] Hopping between different lattices and families
>
> Hi Antonio,
>
> I'm sorry, I don't understand your question. What exactly doesn't work,
> and what do you want to achieve?
>
> Best,
> Anton
>
> On Mon, Oct 31, 2016 at 11:52 AM, <antonio...@usp.br> wrote:
>
>> Hello everyone,
>>
>> I am dealing with exactly the same problem, but in other context (I want
>> to create a one-dimensional p-wave superconductor with two lattices). But
>> the thing is that when I try to implement it just like above, I am
>> neglecting the lat_h(0,y) position, and also losing the lat_e(L,y)
>> position. How can I take this positions without counting two times the
>> "bulk" positions?
>>
>> def make_system(a=1.0, W=1, L=50, mu=1.0, Delta=1.0, t=1.0):
>> # Start with an empty tight-binding system and two square lattices,
>> # corresponding to electron and hole degree of freedom
>> lat_e = kwant.lattice.square(a=1.0, name='e')
>> lat_h = kwant.lattice.square(a=1.0, name='h')
>>
>> sys = kwant.Builder()
>>
>>  Define the scattering region. 
>> sys[(lat_e(x, y) for x in range(L) for y in range(W))] = - mu
>> sys[(lat_h(x, y) for x in range(L) for y in range(W))] = mu
>>
>> # hoppings for both electrons and holes
>> sys[lat_e.neighbors()] = -t
>> sys[lat_h.neighbors()] = t
>>
>> # Superconducting order parameter enters as hopping between
>> # electrons and holes
>> sys[lat_e(x,y), lat_h(x+1,y) for x in range(L-1) for y in range(W)] =
>> Delta
>>
>> Cheers.
>> Antônio.
>>
>> --
>>
>> *De: *"Anton Akhmerov" <anton.akhme...@gmail.com>
>> *Para: *"Camilla Espedal" <camilla.espe...@ntnu.no>
>> *Cc: *kwant-discuss@kwant-project.org, "Bas Nijholt" <
>> basnijh...@gmail.com>
>> *Enviadas: *Segunda-feira, 31 de Outubro de 2016 13:22:01
>> *Assunto: *Re: [Kwant] Hopping between different lattices and families
>>
>>
>> Hi Camilla,
>>
>> It is straightforward to check whether a site is in the system or if its
>> position is inside the shape: http://nbviewer.jupyter
>> .org/url/antonakhmerov.org/misc/site_in_syst.ipynb
>>
>> Basically the site that you want is on the border of the shape and
>> doesn't get includes due to using a strict inequality. As an extra remark I
>> recommend not relying on strict vs non-strict inequality behavior due to
>> float variable rounding. This is especially important if you are dealing
>> with lattices with non-integer unit vectors. Instead a more robust solution
>> would be to have e.g. (-a/2 < x < L+a/2) in the shape definition.
>>
>> Best,
>> Anton
>>
>> On Mon, Oct 31, 2016 at 9:48 AM, Camilla Espedal <camilla.espe...@ntnu.no
>> > wrote:
>>
>>> Thank you!
>>>
>>>
>>>
>>> That worked, but I still don't see why the last thing I did didn't work?
>>>
>>>
>>>
>>> I added a hopping for two points that were within bounds even if
>>> shape_sr started at 1 instead of 0? What am I missing?
>>>
>>>
>>>
>>> Best,
>>>
>>> Camilla
>>>
>>>
>>>

Re: [Kwant] Hopping between different lattices and families

[Anton Akhmerov]

Hi Antonio,

I'm sorry, I don't understand your question. What exactly doesn't work, and
what do you want to achieve?

Best,
Anton

On Mon, Oct 31, 2016 at 11:52 AM, <antonio...@usp.br> wrote:

> Hello everyone,
>
> I am dealing with exactly the same problem, but in other context (I want
> to create a one-dimensional p-wave superconductor with two lattices). But
> the thing is that when I try to implement it just like above, I am
> neglecting the lat_h(0,y) position, and also losing the lat_e(L,y)
> position. How can I take this positions without counting two times the
> "bulk" positions?
>
> def make_system(a=1.0, W=1, L=50, mu=1.0, Delta=1.0, t=1.0):
> # Start with an empty tight-binding system and two square lattices,
> # corresponding to electron and hole degree of freedom
> lat_e = kwant.lattice.square(a=1.0, name='e')
> lat_h = kwant.lattice.square(a=1.0, name='h')
>
> sys = kwant.Builder()
>
>  Define the scattering region. 
> sys[(lat_e(x, y) for x in range(L) for y in range(W))] = - mu
> sys[(lat_h(x, y) for x in range(L) for y in range(W))] = mu
>
> # hoppings for both electrons and holes
> sys[lat_e.neighbors()] = -t
> sys[lat_h.neighbors()] = t
>
> # Superconducting order parameter enters as hopping between
> # electrons and holes
> sys[lat_e(x,y), lat_h(x+1,y) for x in range(L-1) for y in range(W)] =
> Delta
>
> Cheers.
> Antônio.
>
> --
>
> *De: *"Anton Akhmerov" <anton.akhme...@gmail.com>
> *Para: *"Camilla Espedal" <camilla.espe...@ntnu.no>
> *Cc: *kwant-discuss@kwant-project.org, "Bas Nijholt" <basnijh...@gmail.com
> >
> *Enviadas: *Segunda-feira, 31 de Outubro de 2016 13:22:01
> *Assunto: *Re: [Kwant] Hopping between different lattices and families
>
>
> Hi Camilla,
>
> It is straightforward to check whether a site is in the system or if its
> position is inside the shape: http://nbviewer.
> jupyter.org/url/antonakhmerov.org/misc/site_in_syst.ipynb
>
> Basically the site that you want is on the border of the shape and doesn't
> get includes due to using a strict inequality. As an extra remark I
> recommend not relying on strict vs non-strict inequality behavior due to
> float variable rounding. This is especially important if you are dealing
> with lattices with non-integer unit vectors. Instead a more robust solution
> would be to have e.g. (-a/2 < x < L+a/2) in the shape definition.
>
> Best,
> Anton
>
> On Mon, Oct 31, 2016 at 9:48 AM, Camilla Espedal <camilla.espe...@ntnu.no>
> wrote:
>
>> Thank you!
>>
>>
>>
>> That worked, but I still don't see why the last thing I did didn't work?
>>
>>
>>
>> I added a hopping for two points that were within bounds even if shape_sr
>> started at 1 instead of 0? What am I missing?
>>
>>
>>
>> Best,
>>
>> Camilla
>>
>>
>>
>>
>>
>> *From:* Bas Nijholt [mailto:basnijh...@gmail.com]
>> *Sent:* 31. oktober 2016 14:41
>> *To:* Camilla Espedal <camilla.espe...@ntnu.no>; Anton Akhmerov <
>> anton.akhme...@gmail.com>
>>
>> *Cc:* kwant-discuss@kwant-project.org
>> *Subject:* Re: [Kwant] Hopping between different lattices and families
>>
>>
>>
>> The error message says it all :-)
>>
>>
>>
>> Using `return ((0 <= x < L) and (0 <= y < W)) ` should fix your problem.
>>
>>
>>
>> Best, Bas
>>
>>
>>
>> On Mon, 31 Oct 2016 at 14:29 Camilla Espedal > <camilla%20espedal%20%3ccamilla.espe...@ntnu.no%3e>> wrote:
>>
>> Hi again,
>>
>> Thanks, I will try to write it more clearly. So the script I use to make
>> the system is:
>>
>> L = 100
>> W = 40
>> a = 1
>> t = 1
>>
>> def make_system(W,L,a,t):
>>
>> # Make the lattices
>> # We define two lattices (up and down) with two sublattices A and B
>> lat_up = kwant.lattice.general([(a,a),(a,-a)],[(0,0),(a,0)], name='up')
>> A_up, B_up = lat_up.sublattices
>> lat_down = kwant.lattice.general([(a,a),(a,-a)],[(0,0),(a,0)],
>> name='down')
>> A_down, B_down = lat_down.sublattices
>>
>> # Define the shape of the scattering region. Must return true where there
>> are sites.
>> def shape_sr(pos):
>> x, y = pos
>> return ((0 < x < L) and (0 < y < W))
>>
>>
>> sys = kwant.Builder()
>>
>> sys[lat_up.shape(shape_sr, (1,1))] = 2
>> sys[lat_down.shape(shape_sr, (1,1))] = 2
>>
>> sys

Re: [Kwant] The tight-binding Hamiltonian in real space and momentum space

[Anton Akhmerov]

Hi Johnny,

To obtain a tight-binding Hamiltonian from a momentum space
Hamitlonian we need to calculate a tight-binding description of the
momentum operator. This is done using discretization, as shown e.g.
here: 
https://kwant-project.org/doc/1/tutorial/tutorial1#discretization-of-a-schrodinger-hamiltonian
or here: 
http://mybinder.org/repo/kwant-project/kwant-tutorial-2016/notebooks/2.1.Discretizing-Hamiltonians-computing-spectra.ipynb

Best,
Anton

On Tue, Nov 1, 2016 at 9:45 AM, T.C. Wu  wrote:
> Dear Kwant users,
>
> I recently reading the course materials about topological insulator at
> http://topocondmat.org/w7_defects/ti_majoranas.html. However, the Kwant code
> which simulates the majoranas on the quantum spin-Hall edge is not quite
> transparent to me. In particular, the hoppings defined in Kwant are
>
> def onsite(site, p):
> (x, y) = site.pos
> return (p.M - 4 * p.B) * s0szsz - 4 * p.D * s0s0sz + p.gaps(x, y)[1]
> * mysz + p.gaps(x, y)[0] * s0s0sx
>
> def hopx(site1, site2, p):
> return p.B * s0szsz + p.D * s0s0sz + 0.5j * p.A * szsxsz
>
> def hopy(site1, site2, p):
> return p.B * s0szsz + p.D * s0s0sz - 0.5j * p.A * s0sysz
>
> I got a bit lost since I don't quite understand why should hopx and hopy be
> defined in this way. Maybe a silly question, but how should one obtain the
> hoppings by looking at the momentum space BdG Hamiltonian?
>
> Thanks for your help.
>
> Best,
> Johnny

Re: [Kwant] Hopping between different lattices and families

[Anton Akhmerov]

Hi Camilla,

Please double-check the error message that you see.

Your assumption why the code doesn't work is not right: it's possible
to add a hopping from any site to any site, regardless of distance or
lattices involved. My best guess is that the sites aren't present in
the system yet.

As a general advice, when describing a problem try to provide complete
information required to reproduce this problem. A script and the error
message would be usually useful.

Best,
Anton

On Mon, Oct 31, 2016 at 6:58 AM, Camilla Espedal
<camilla.espe...@ntnu.no> wrote:
> Hi again,
>
> I tried to add
>
> sys[A_up(2,2), B_down(2,2)] = 2
>
> This does not work, and I think it is because A_up and B_down are not only on 
> different sublattices, but on different lattices as well. In the tutorial on 
> superconductors 
> (https://kwant-project.org/doc/1.0/tutorial/tutorial5#lattice-description-using-different-lattices)
>  they define the hopping as
>
> sys[((lat_e(x, y), lat_h(x, y)) for x in range(Deltapos, L)
>  for y in range(W))] = Delta
>
> between the lattices but on the same spatial point. But I want to hop between 
> two different lattices from one point in space to another, if that makes 
> sense.
>
> Best,
> Camilla
>
> -Original Message-
> From: Anton Akhmerov [mailto:anton.akhme...@gmail.com]
> Sent: 31. oktober 2016 11:48
> To: Camilla Espedal <camilla.espe...@ntnu.no>
> Cc: kwant-discuss@kwant-project.org
> Subject: Re: [Kwant] Hopping between different lattices and families
>
> Hi Camilla,
>
> It's exactly like you would expect: syst[A_up(i, j), B_down(i, j)] = value. 
> See e.g. https://kwant-project.org/doc/1/tutorial/tutorial4
>
> Best,
> Anton
>
> On Mon, Oct 31, 2016 at 6:43 AM, Camilla Espedal <camilla.espe...@ntnu.no> 
> wrote:
>> Hi,
>>
>>
>>
>> To explain what I mean. I have a system where I have separated spin-up
>> and spin-down into two lattices (like electron and hole in the example
>> at the kwant site) so that it is would be easier to extract
>> spin-resolved information (G_up/up, Gup/down etc.). In addition, these
>> two lattices consists of two sublattices (A and B).  The code is like
>> this
>>
>>
>>
>> lat_up =
>> kwant.lattice.general([(1,0),(s,c)],[(0,0),(0,1/sqrt(3))],
>> name='up')
>>
>> A_up, B_up = lat_up.sublattices
>>
>> lat_down =
>> kwant.lattice.general([(1,0),(s,c)],[(0,0),(0,1/sqrt(3))],
>> name='down')
>>
>> A_down, B_down = lat_down.sublattices
>>
>>
>>
>> What I want to do now, is to add a hopping term between say A_up (i,j)
>> and B_down (i.j). So a hopping term between the two different lattices
>> and sublattices. How can I implement this? Or is there a better way to
>> achieve what I want?
>>
>>
>>
>> Best
>>
>> Camilla Espedal

Re: [Kwant] Hopping between different lattices and families

[Anton Akhmerov]

Hi Camilla,

It's exactly like you would expect: syst[A_up(i, j), B_down(i, j)] =
value. See e.g. https://kwant-project.org/doc/1/tutorial/tutorial4

Best,
Anton

On Mon, Oct 31, 2016 at 6:43 AM, Camilla Espedal
 wrote:
> Hi,
>
>
>
> To explain what I mean. I have a system where I have separated spin-up and
> spin-down into two lattices (like electron and hole in the example at the
> kwant site) so that it is would be easier to extract spin-resolved
> information (G_up/up, Gup/down etc.). In addition, these two lattices
> consists of two sublattices (A and B).  The code is like this
>
>
>
> lat_up = kwant.lattice.general([(1,0),(s,c)],[(0,0),(0,1/sqrt(3))],
> name='up')
>
> A_up, B_up = lat_up.sublattices
>
> lat_down = kwant.lattice.general([(1,0),(s,c)],[(0,0),(0,1/sqrt(3))],
> name='down')
>
> A_down, B_down = lat_down.sublattices
>
>
>
> What I want to do now, is to add a hopping term between say A_up (i,j)  and
> B_down (i.j). So a hopping term between the two different lattices and
> sublattices. How can I implement this? Or is there a better way to achieve
> what I want?
>
>
>
> Best
>
> Camilla Espedal

Re: [Kwant] kkghosh

[Anton Akhmerov]

Dear Kamal,

It appears that your questions mostly relate to the background physics
knowledge and are not specific to Kwant. There is vast literature
explaining how quantum transport works, and it provides a much more
detailed and systematic answer to your questions than this mailing
list can do.

Best regards,
Anton Akhmerov

On Fri, Oct 28, 2016 at 12:44 PM, kamal ghosh <kk_gh...@rediffmail.com> wrote:
> Hi all,
> Many thanks to all the developers of kwant not only for the strong potential
> development of the code but
> very specially also for the help to solve coding and related physics
> problems. Presently, I am wondering
> how to calculate current in a device, say a QWR, connected to source and
> drain through two infinite
> leads! Applying the Landauer relation, I require no. of channels M upto the
> level of Fermi and
> Transmission probabilities per channel. To me, this appears very involved
> programming; because (i) first
> to find the Fermi first and next to find the no. of channels below that (ii)
> secondly, for any channel
> the all energies upto the Fermi I have to calculate and (iii) finally to sum
> up. Then we may get the
> total transmission probability.
> I am not sure that is it right to think of a constant probability in a
> particular channel? I think it is
> "no".
> However, I noticed a tutoring correspondence from Joe and Akhmerov in a
> correspondence with an user
> working with calculating the current like me. I could not understand the
> solution because current carried
> by scattering state is written by a code which I am afraid how to write that
> code in high level python
> text? Please, bear with me on this trifling matter. Should we wait for the
> version kwant2?
> All the best wishes and awaiting your wise advise.
> K.K.Ghosh

Re: [Kwant] barrier transport

[Anton Akhmerov]

Dear Kamal,

* There already is a tutorial showing transmission through a barrier:
https://kwant-project.org/doc/1/tutorial/tutorial2#spatially-dependent-values-through-functions
* Looking at your script I see that you define the function
"plot_transmission" but don't actually call it. Could this be the
origin of your problem?

Best,
Anton

On Sun, Oct 23, 2016 at 7:15 PM, kamal ghosh  wrote:
> Dear All,
>
> I have tried with a kwant code to see the transmission through barrier
> imposed by a qwr between the
> source and drain. The wire is constructed and its bands are formed
> successfully; but the plot of
> transmission through the constant barrier is not obtained. In the long run,
> I would like to use the
> standardised varying barrier potential and see the transport.
> Please help me overcome this difficulty so that I may proceed further.
> Kind regards.
> K.K.Ghosh
> PS. code is attached herewith.
>

Re: [Kwant] characteristic graphs

[Anton Akhmerov]

Dear Kamal,

> (i) as to how to draw the VD- ID, VG-ID etc. by kwant coding. I am to some
> extent clear to
> get those relations from the back calculations of G= ID/VG and the current
> expression in
> terms of the G. Finding G by kwant they may be found out. But, I like very
> much a direct
> coding technique to get these curves.

I believe that getting a differential conductance using Kwant and then
integrating it numerically is the only approach in this case.

> (ii) as to how to visualize a 3D graphics of nanostructure like QWR, QD
> which I visualize in
> 2D without problem.

You can either use kwant.plot using a 3D system, or any other 3D
plotting library if the capacity of kwant.plot is insufficient.

Best regards,
Anton Akhmerov

Re: [Kwant] Runnning Kwant computations on parallel machines

[Anton Akhmerov]

> The common approach to simulating mesoscopic transport in large
> devices is to use an effective tight-binding model instead of physical
> atoms. That way you can efficiently simulate large devices still on
> your laptop. While the approach is well known and broadly used since
> the beginning of graphene days, there were two recent publications
> dedicated specifically to making such a simulation match an experiment
> in graphene devices [1]. I believe that approach is the optimal way
> for you to proceed.

Sorry, forgot to add the references; here they go:
http://arxiv.org/abs/1605.06924 and http://arxiv.org/abs/1407.5620


> Best,
> Anton
>
> On Wed, Aug 31, 2016 at 10:50 AM, Christoph Groth
>  wrote:
>> Adrian Nosek wrote:
>>
>>> My name is Adrian Nosek, I am a graduate student at the University of
>>> California Riverside in Condensed Matter Physics, and I was wondering
>>> whether there are any plans to extend Kwant and let it run on parallel
>>> machines, so that systems of bigger sample size could be modeled?
>>
>>
>> Hi Adrian, thanks for your interest.  Parallel computations are usually
>> subdivided into two broad classes that can be combined:
>>
>> • Task parallelism: Spread independent calculations (e.g. at  different
>> energies, or for different disorder realizations)  across different
>> cores/machines.
>>
>> • Data parallelism: Spread a single computation.  For example,  each
>> core/machine could take care of part of a huge Hamiltonian  matrix.
>>
>> There are many variants of both of course.  Each can be done across a single
>> machine with many cores, across an interconnected cluster of machines, or
>> across machines in a "grid" or the "cloud".
>>
>> Task parallelism can be done with Kwant in various ways.  One way that we
>> recommend (concurrent.futures) is illustrated with an complete example (the
>> spin valve) in Kwant paper.
>>
>> Basic data parallelism can be achieved by compiling Kwant with a parallel
>> BLAS/LAPACK library.  The default Debian/Ubuntu packages of Kwant, for
>> example, are built with the parallel OpenBLAS. (Unfortunately, with OpenBLAS
>> the speedup seems to be negligible.) But perhaps it’s a matter of tuning
>> OpenBLAS properly.  It’s also possible that another BLAS (e.g. intel MKL)
>> would work better.
>>
>> There is so far no Kwant solver that utilizes data parallelism other than
>> through BLAS, but it should be not too difficult to adapt the (currently
>> sequential) MUMPS solver to use MUMPS with MPI.  There have been some
>> discussions about this on the kwant-devel list:
>>
>> https://www.mail-archive.com/kwant-devel%40kwant-project.org/msg00047.html
>> https://www.mail-archive.com/kwant-devel%40kwant-project.org/msg00052.html
>>
>> We would be of course happy to support anyone who would like to work on
>> this!
>>
>>> I attempted to calculate the conductivity of a sample with 20 atoms (I
>>> believe that's what it was but its been a while since I ran the simulation)
>>> and this is where my machine ran out of RAM.
>>
>>
>> How many orbitals per site does your model feature?  How much RAM does your
>> machine have?  On my laptop with 8 GB of RAM, I’ve done calculations with up
>> to 2 millions of single-orbital sites (see the benchmark in the Kwant
>> paper).
>>
>> But it’s true that the MUMPS solver that Kwant uses by default (and which is
>> the fastest) needs a lot of RAM.
>>
>>> I really like the idea of simulating my experiment using Kwant, but is
>>> there any advantage to model let's say a 4 by 20 micron size system using
>>> Kwant vs a classical approach? I believe the mean free paths of electrons in
>>> graphene can easily reach micron-long pathways in high mobility graphene-hBN
>>> sandwiches, so that is the reason why I am asking.
>>
>>
>> Coherent simulations (like the ones done by Kwant) make sense when the
>> coherence length is comparable to the system size.  Whether it makes sense
>> for your particular experiment I don’t know.
>>
>>> Also, would it make any sense to implement such a simulation in the
>>> context of Big Data Analysis, where analysis is performed on parallel
>>> working machines where one could diagonalize such a big Hamiltonian, e.g.
>>> using Apache Spark?
>>
>>
>> Can Apache Spark be used for solving huge systems of linear equations?
>> Because this is what Kwant is doing internally.  My impression is that "big
>> data" frameworks are rather meant for problems that require less
>> communication.
>>
>> But there are certainly libraries for solving linear systems in parallel,
>> see the threads linked above.