Re: [Kwant] [EXTERNAL] An issue with plotting two systems with different dimension in one system

[Joseph Weston (Aquent LLC - Canada)]

Hi,

The problem is that your ‘lat’ lattice exists in R^3, whereas your ‘honeycomb’ 
lattice exists in R^2. Kwant won’t try to guess what embedding you want to use 
for the ‘honeycomb’ lattice, so if you want to use both in the same system you 
will have to specify your honeycomb lattice using vectors in R^3.

If, for example, you wanted to embed your honeycomb lattice in the x-y plane 
you would do:

honeycomb = kwant.lattice.general([(1, 0, 0), (1 / 2, sqrt(3) / 2, 0)], [(0, 0, 
0), (0, 1 / sqrt(3), 0)])

Note the extra zeros compared to the vectors provided in your code.

It’s true that the error message is not very informative in this case. There is 
actually an issue already open about this: 
https://gitlab.kwant-project.org/kwant/kwant/-/issues/222, however nobody yet 
took the time to fix it.

Happy Kwanting,

Joe

From: Kwant-discuss  On Behalf Of 
tavakkolidjawad
Sent: Sunday, May 3, 2020 9:10 AM
To: kwant-discuss@kwant-project.org
Subject: [EXTERNAL] [Kwant] An issue with plotting two systems with different 
dimension in one system


dear all

I want to plot two different structures in one system. Although it is easy to 
plot separate lattices with the same structure, I will not be able to solve the 
problem when the structures are different. For example, I plotted two lattices 
with the same diamond structure, but when I replaced one of them with a 
honeycomb structure, I got an error after plotting the system. Apparently and 
logically, I think there is no bug in my code.

Is there a way to solve this problem?

I ask you to take a look at my code and report any possible bugs to me.

Thanks



###

import kwant
from math import sqrt

lat = kwant.lattice.general([(0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0)],
[(15, 15, 15), (15.25, 15.25, 15.25)], name=['a1', 'a2'])
a1, a2 = lat.sublattices



honeycomb = kwant.lattice.general([(1, 0), (1 / 2, sqrt(3) / 2)],
[(0, 0), (0, 1 / sqrt(3))])
subA , subB = honeycomb.sublattices



def make_system(a=25, b=25, c=30):

syst = kwant.Builder()

def cuboid_shape(pos):
x, y, z = pos
return 15 <= x < a and 15 <= y < b and 15 <= z < c

syst[a1.shape(cuboid_shape, (15, 15, 15))] = 1
syst[a2.shape(cuboid_shape, (15, 15, 15))] = 1
syst[lat.neighbors()] = 1


def graphene_shape(pos):
x,y = pos
return 0 <= x < 10 and 0 <= y < 10

syst[honeycomb.shape(graphene_shape,(0,0))] = 1
syst[honeycomb.neighbors()] = 1


return syst



def main():

syst= make_system()

kwant.plot(syst)


syst = make_system(a=1.1, b=1.1, c=1.1)

if __name__ == '__main__':
main()

###

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

[Joseph Weston (Aquent LLC - Canada)]

Hi Anton,

Thanks for the detailed response.

I think I understand most of your answer, but I'm a little confused about one 
part.

If I remember correctly, Kwant computes self-energies for leads using the mode 
decomposition via V†ΦΛΦ^-1, where V is the hopping, Φ Is the matrix of outgoing 
modes (Φ^-1 its right-inverse) and Λ is the matrix of translation eigenvalues.
The only places I can see where a pole may appear infinite translation 
eigenvalues λ or badly conditioned Φ.
On the other hand, when constructing the linear system for the scattering 
problem we extend the left-hand-side with terms like V†ΦΛ.

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.

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 amplitudes, and 
> the RHSs are indicator vectors with the 1’s in this “extended” part, which 
> corresponds to a “single incoming mode” bou

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

[Joseph Weston (Aquent LLC - Canada)]

Hello again,

I Noticed a couple of minor mistakes in my previous email:


  1.  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.
  2.  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 amplitudes, and the 
RHSs are indicator vectors with the 1's in this "extended" part, which 
corresponds to a "single incoming mode" boundary condition.

It seems to me that the salient difference is in the boundary conditions, and I 
do not have a good intuition as to why one set of boundary conditions would 
make the linear system easier/harder to solve.

Happy Kwanting!

Joe

[1]: 
https://gitlab.kwant-project.org/kwant/kwant/-/blob/master/kwant/solvers/common.py#L428<https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgitlab.kwant-project.org%2Fkwant%2Fkwant%2F-%2Fblob%2Fmaster%2Fkwant%2Fsolvers%2Fcommon.py%23L428=02%7C01%7Cv-josewe%40microsoft.com%7Cccd4df85ea5a44fb98e508d7cb6fbfa3%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637201553194781503=jNV1lZwk4xs%2BfX25KOuP0Z992wS3as9Uribwou1P1rk%3D=0>
[2]: 
https://iopscience.iop.org/article/10.1088/1367-2630/16/6/063065/pdf<https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fiopscience.iop.org%2Farticle%2F10.1088%2F1367-2630%2F16%2F6%2F063065%2Fpdf=02%7C01%7Cv-josewe%40microsoft.com%7Cccd4df85ea5a44fb98e508d7cb6fbfa3%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637201553194791491=5ltAdppLxjv%2BAf6fHEb4xTlSrZ%2FcGNY1jPFwnWO932s%3D=0>
[3]: 
https://epub.uni-regensburg.de/12142/<https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fepub.uni-regensburg.de%2F12142%2F=02%7C01%7Cv-josewe%40microsoft.com%7Cccd4df85ea5a44fb98e508d7cb6fbfa3%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637201553194791491=B1HTfGTAeQqmX%2BuLFsgxrGdtMYMij%2BbDkgupLe6efOg%3D=0>

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

[Joseph Weston (Aquent LLC - Canada)]

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 amplitudes, and the 
RHSs are indicator vectors with the 1's in this "extended" part, which 
corresponds to a "single incoming mode" boundary condition.

It seems to me that the salient difference is in the boundary conditions, and I 
do not have a good intuition as to why one set of boundary conditions would 
make the linear system easier/harder to solve.

Happy Kwanting!

Joe

[1]: 
https://gitlab.kwant-project.org/kwant/kwant/-/blob/master/kwant/solvers/common.py#L428
[2]: https://iopscience.iop.org/article/10.1088/1367-2630/16/6/063065/pdf
[3]: https://epub.uni-regensburg.de/12142/

Re: [Kwant] Problems to add disorder to the Kitaev chain

[Joseph Weston]

Hi Luca,


> thanks again for your useful reply.
>
> If I understand correctly, then I am not sure
> how to simulate the physics of disorder.
>
> Indeed, that task requires to create various
> (in principle infinite in number) configurations,
> measure some physical quantity on each of them, and then
> to take the average of all these measurements on the created configurations.
>
> How can I create the different configurations ?
> Starting e.g. from the piece of code

Indeed, typically you'll want to take the average over different
disorder configurations. This means using different values for the
'salt' parameter. The actual values used for 'salt' are not important,
as long as they differ; 'kwant.digest.uniform' is a hash function (see
https://en.wikipedia.org/wiki/Salt_(cryptography) ).

You would do something like the following:


    syst = kitaev(L=...).finalized()

    params = dict(mu=..., W=..., t=..., Delta=...)

   

    # We (arbitrarily) choose the salts as the strings '0', '1', ...

    salts = map(str, range(100))

    for params['salt'] in salts:

        ham = syst.hamiltonian_submatrix(params=params, sparse=False)

        evals = np.linalg.eigvalsh(ham)

        ...


Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] What is the algorithm for calculating the mode decomposition

[Joseph Weston]

Me again,


I just noticed that the notebook at http://tiny.cc/kwant-journal-club
looks kind of wonky because nbviewer doesn't understand some of the
markup I used. I've now made the notebook available on binder as well
(even though there's no executable code really):

http://tiny.cc/kwant-journal-club-notes

Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] What is the algorithm for calculating the mode decomposition

[Joseph Weston]

Hi Zhan,


>
> In the latest Kwant manual describing the details about
> kwant.physics.modes, I find a sentence "This function uses the most
> stable and efficient algorithm for calculating the mode decomposition
> that the Kwant authors are aware about. Its details are to be
> published." Do anyone know what the algorithm is? I look forward some
> references if possible.


That is an excellent question.

There is a review of numerical quantum transport currently in
preparation by the main Kwant authors, which will include explicit
details of the algorithm implemented by Kwant. Unfortunately the review
is not yet in a state where it can be shared, but rest assured that as
soon as it is we will inform the Kwant community via the mailing list.


In the meantime there are a few places where you can get some more
information:


1. https://tiny.cc/maryland-kwant-tutorial

   This is a tutorial I gave at the university of Maryland about how
Kwant works. It gives a high-level overview of the scattering problem
and the lead-modes problem, but doesn't go into details though.


2. http://tiny.cc/kwant-journal-club

    This is a talk I gave at Delft university that dives deep into the
details. It contains code snippets from the Kwant source code, and links
back to the relevant places in the source code itself.


3. https://gitlab.kwant-project.org/kwant/kwant/

    The Kwant source code itself! This is of course the canonical place
to find out how Kwant works internally. The start of the modes-related
code is in 'kwant/physics/leads.py' around line 990, but it is quite
dense (and that's putting it mildly!).


Hopefully that will be enough to get you started; let us know how you
get on.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Calculating the total energy in a superconducting lead

[Joseph Weston]

Hi Jannis,


> It does make sense and helps a lot.
>
> Is there a way to calculate E_n(k) of a 2D infinite system in kwant?
> 2D translational symmetric systems can't be finalized and
> cell_hamiltionian and inter_cell_hopping don't work with wraparound
> either.
>

This is something that we've wanted for a while, but now it's on our
immediate roadmap (see [1]).


Wraparound returns a Builder with *no* translational symmetries, but
with extra momentum parameters, so you should use
'hamiltonian_submatrix' in this case.

'syst.hamiltonian_submatrix(params=dict(k_x=0.5, k_y=0.6))' directly
gives you the k-space Hamiltonian H(k_x, k_y) at the requested momentum
when 'syst' is a system finalized from a Builder returned from wraparound.


Happy Kwanting,

Joe


[1]: https://gitlab.kwant-project.org/kwant/kwant/issues/316




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] matrix

[Joseph Weston]

Hi Nafise,


> Dear all,
>
> Is there any one to know how we can use one element of matrix? For
> example we have matrix A as follows:
>
> A=[2 4;6,0]
>
> We have A matrix with dimension 2*2. If  we want to use one element of
> matrix such as A(2,1), (that is number 6 in the matrix), How we can
> write in kwant? Any help appreciate.


Well it depends what you mean.

In Kwant you typically construct your system using a Builder. You can
think of this as *analogous* to a sparse matrix where you can index the
entries with *pairs of sites* (hoppings) instead of integers.

After finalizing a builder you can call 'hamiltonian_submatrix', which
will give you an *actual* dense or sparse matrix (depending on whether
you pass 'sparse=True' or 'sparse=False'). You can then index this
matrix as you would normally.


Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Calculating the total energy in a superconducting lead

[Joseph Weston]

Hi Jannis,


>
> thank you for the fast response. You are probably right, that I need
> to use the k-space Hamiltonian for the infinite case. Would my method
> work for a finite system then? Potentially with periodic boundary
> conditions?


For a finite system (with or without PBC) you should directly
diagonalize Hamiltonian and count the number of states below the Fermi
energy.


> The measure, that I am looking for, is the total energy of all
> occupied particle states. I am not really sure, how to calculate this
> in an infinite system. First I probably need to isolate the particle
> Hamiltonian with a projector. But I don't know how to proceed from
> there in k-space.

If I understand your question correctly then you need to do do a k-space
integration; something like:


 ∑_n ∫ E_n(k) f(E_n(k)) dk


Where the sum runs over all the bands, the integral runs over the
Brillouin zone and 'f' is the occupation at energy E (e.g. Fermi-Dirac).
In this way you count the energy contribution for each state, weighted
by its occupation. Note that we have elided the density of states, as we
are integrating directly in momentum space, as opposed to in energy, and
the density of states in 1D is constant in k.

Does that make sense?


Aside from this I have the suspicion that you might not necessarily get
the result you are looking for. You say that you are dealing with an
s-wave superconductor, but want to calculate the total energy of the
electrons only. However in an s-wave superconductor the eigenstates are
not electrons and holes, but a *superposition* of electron and hole (the
so-called Bogoliubov quasiparticles). I would therefore recommend taking
care to look at what it is that you want to calculate, and asking
whether it is meaningful.


Happy Kwanting,

Joe

∫ {\displaystyle \displaystyle \int } 



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Calculating the total energy in a superconducting lead

[Joseph Weston]

Hi again,

I noticed my code example probably won't even run in the generic case;
we need to do this wrangling with 'numpy.zeros_like' because
'inter_cell_hopping' returns a rectangular matrix in general (when only
a subset of the sites in a unit cell connect via hoppings to the next
unit cell). Probably the code should look something like this:


        def H_k(k, params):

        H_0 = syst.cell_hamiltonian(args, params=params)
   
    hop = syst.inter_cell_hopping(args, params=params)
        V = np.zeros_like(H_0)
    V[:, : hop.shape[1]] = hop

    return H_0 + V * exp(-1j * k) + V.conjugate().transpose() *
exp(1j * k)


Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Calculating the total energy in a superconducting lead

[Joseph Weston]

Hi Jannis,


You mentioned that you system is infinite and has translational
symmetry, yet you are using 'hamiltonian_submatrix'.

'hamiltonian_submatrix' probably does not return what you want; what
should it return when the system is infinite, as it is in your case?
Probably you want to construct the k-space Hamiltonian and work with
that. You could do that with the following function:


    def H_k(k, params):

    H_0 = syst.cell_hamiltonian(params=params)

    V = np.zeros_like(H_0)

    V += syst.inter_cell_hopping(params=params)

    return H_0 + V * exp(1j * k) + V.conjugate().transpose() *
exp(-1j * k)


If you are just interested in the eigenvalues/vectors you can also use
'kwant.physics.Bands', which will allow you to evaluate the
eigenvalues/vectors given a k (check out the documentation to find out
how to get the eigenvectors in addition to the eigenvalues).

I would first rectify this before looking any further.

Happy Kwanting,

Joe


> I am trying to calculate the total energy of the electrons in an
> s-wave superconducting lead with spins. But I am not sure, if my
> method is correct.
>
> Currently I separate the particle and hole Hamiltonians with a projector
>
>     H = syst.hamiltonian_submatrix(params=params)
>     P_electron = np.kron(np.eye(L+1), np.array([[1, 0, 0, 0], [0, 1,
> 0, 0], [0, 0,0,0], [0, 0,0,0]]))
>     H_electron = P_electron.conjugate().transpose() @ H @ P_electron
>
> where L is the number of sites. L+1 because I get a dimension mismatch
> otherwise. My guess here is, that kwant generates an extra site
> because of the translational symmetry. For the conservation laws I used
>
>     particle_hole=pauli.sxs0, conservation_law=-pauli.szs0
>
> where pauli.sisj is the kroneker product of the i-th and the j-th
> pauli matrix. For my Hamiltonian I used the base (e_up, e_down,
> h_down, h_up).
>
> Then I calculate the eigenvalues of H_electron and multiply them with
> the fermi distribution to get the total energy.
>
>     spectrum=  np.linalg.eigvalsh(H_electron)
>     fermi_dist=kwant.kpm.fermi_distribution(spectrum, params["mu"], T)
>     total_energy=np.dot(fermi_dist,spectrum)
>
> Now to the part, why I suspect a mistake in my method. I use a
> spin-dependent hopping term and scan over a parameter /b/, that
> changes the spin dependency. When doing so the qualitative behavior of
> the total energy /E_total(b)/ changes with the size of my unit cell.
> But all sites in my unit cell are identical. So I would suspect, that
> the results should only differ by a factor. For large systems the
> curvature of /E_total(b)/ converges, but for small systems the results
> vary quit a lot.
>
> I'd really appreciate some help on this.
> Best regards,
> Jannis
>


signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Access to eigenvalue, eigenvector and number of points in each unit cell

[Joseph Weston]

Hi Sajad,


> Dear all,
>
> I need to access to the number of atoms of my unit cell,
> the eigenvalue and eigenvectors for each eigenvalue  of my Hamiltonian.
>
> Is there any one to help me and let me know if it is possible in kwan
> to access them.
>

Could you post a short code example showing what you are doing? You
refer to "atoms" and "unit cell", so I presume that you are trying to
create a system with translational symmetry, and that each Kwant *site*
corresponds to a single atom.


It is not 100% clear to me what you want when you say "the eigenvectors
and eigenvalue" of your Hamiltonian; if your system has translational
symmetry then presumably you want the eigen-decomposition *at a given
quasi-momentum*, but you do not explicitly state this, so I am not sure.


Posting a complete code example is useful because it is more precise
than describing your problem with words.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Non-equilibrium Green's functions and Current in tkwant

[Joseph Weston]

Hi Rudolf,


> The kwant documentation on kwant.solvers.common.GreensFunction() states
> that calling this method will return the retarded Green's function. Does
> this method evaluate the equilibrium or non-equilibrium retarded Green's
> function? I am interested in the time-dependent non-equilibrium density
> matrix, which can be written as the kinetic lesser NE Green's function
> G_<(t,t')
>
> rho_neq (t) = -iG_<(t,t')|t=t'.


'kwant.greens_function' essentially gives you a few matrix elements of
the following quantity: Gᴿ(E) = (H + Σ(E))⁻¹ where Σ is the self-energy
of the leads at energy E. Specifically it gives you the matrix elements
that connect the interface sites of the different leads; these are the
matrix elements needed to calculate transport quantities (like conductance).

In order to calculate the lesser Green's function for the
time-independent, nonequilibrium case (i.e. different chemical
potentials / temperatures in different leads) you would need to take a
different approach I think.

One way to do so would be to apply eq. 22 of
https://arxiv.org/abs/1307.6419 to the special case where Ψ_αE(t) has
trivial time dependence: Ψ_αE(t) = exp(-iEt) Ψ_αE where Ψ_αE is the
scattering state corresponding to incoming mode α at energy E (i.e. what
kwant.wave_function gives you).

If your system is *time-dependent* then you will need to calculate the
time-evolved scattering states, which Kwant cannot give you directly.

Indeed I did some work on calculating such quantities during my PhD.,
and one of the outputs was a code "tkwant", which is essentially a
solver for  working with Kwant systems that include time-dependence. It
essentially does the time evolution and energy integration for you.

tkwant is still actively developed, but I am no longer directly
involved; you can check it out here:
https://gitlab.kwant-project.org/kwant/tkwant/ there's also links to the
documentation which includes some examples which should help you get
started if this is the direction you want to go.


Happy (t)kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Josephson Junctions using tkwant

[Joseph Weston]

Hi again Denise,

I'm just writing again to clarify something I wrote in my previous reply:

> Unfortunately tkwant does not presently have a maintainer, and no
> developer time has been allocated to it. This may change in the future,
> but at this point we cannot provide any promises. tkwant was anyway an
> experimental code that is not "production ready" in the same way that
> kwant is.


This is not actually true; there has recently been quite a bit of work
on tkwant recently, and Thomas Kloss is the current maintainer. I was
aware of this, but have been out of the loop on tkwant development,
which lead to me writing the misleading statement quoted above.

If the tutorials do not work with the current version of tkwant then
this is a bug (either in the documentation or tkwant itself), and I
would encourage you to open an issue on the tkwant issue tracker [1].

Sorry again for the misleading message.

Happy tkwanting,

Joe


[1]: https://gitlab.kwant-project.org/kwant/tkwant/issues/new




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Josephson Junctions using tkwant

[Joseph Weston]

Hi Denise,

>
> Has anyone recently tried to simulate Josephson Junctions using tkwant?
> I tried to follow the examples and tutorials, but several of them do
> not work in the last tkwant version.


Unfortunately tkwant does not presently have a maintainer, and no
developer time has been allocated to it. This may change in the future,
but at this point we cannot provide any promises. tkwant was anyway an
experimental code that is not "production ready" in the same way that
kwant is.


> I am sending attached and you can find below the code that I am
> running but the current does not oscillate as a function of time and
> the IV curve is linear. Any suggestions on what may be wrong? Any help
> or comment is appreciated.

Just skimming your code I notice the following:


>    
>     ### add JJ phase
>     syst[(lat(-1,0),lat(0,0))] =  -t*np.exp(1j*phase(time))*c_e
> -t*np.exp(-1j*phase(time))*c_h  ### JJ phase
>    


Here you are assigning a *constant value* to the hopping. Presumably you
would like instead to have a *function* that depends on time. You should
instead write something like:


    def hopping(site_a, site_b, time):

        return -t*np.exp(1j*phase(time))*c_e -t*np.exp(-1j*phase(time))*c_h


    syst[(lat(-1,0),lat(0,0))] = hopping


Happy (t)kwanting,


Joe





signature.asc
Description: OpenPGP digital signature

Re: [Kwant] question

[Joseph Weston]

Hi Nafise

>
> I need to make a periodic lattice with hole. In fact I should make holes on 
> the scattering region and also on the leads. Although I can make this kind of 
> lattice by kwant, I have problem about the distances between holes. I want to 
> make a periodic holes on the nanoribbon but the distance between holes in the 
> scattering region is different from the distance between holes in the other 
> regions. Would you please let me know How Can I make a same distance between 
> holes? Should I work on the translational symmetry in the leads?

Because the leads need to be translationally invariant if you want a
very large distance between defects then your unit cell in the leads
needs to be correspondingly large.

If you post a code example of what your problem is specifically we may
be able to help more.


Happy Kwanting,

Joe


P.S. sorry for the double reply; I forgot to send to the mailing list also




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] help conda

[Joseph Weston]

Hi,

> Dear Sir,
> I have a problem to install Kwant when download the Pre-built packages ( 
> Anaconda package ). That is to say I can't download the package. My  
> operating systems is Microsoft Windows. Can you help me.


Please provide more information; what exact steps did you take when
trying to install Kwant. Did you use the command line and the following
command:

conda install -c conda-forge kwant

or did you use the Anaconda GUI to install the package? In the latter case did 
you make sure to enable the "conda-forge" channel?

Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] How to solve the Hamiltonian matrix with off-diagonal terms?

[Joseph Weston]

Hi,

> We consider the values of the  conservation law based on the
> off-diagonal term of Hamiltonian matrix, details are as follows:
>
>  
>
> H = tinyarray.array([[0, 0, 0, 1, 0, 0],
>
>     [0, 0, 0, 0, 1, 0],
>
>     [0, 0, 0, 0, 0, 1],
>
>     [1, 0, 0, 0, 0, 0],
>
>     [0, 1, 0, 0, 0, 0],
>
>     [0, 0, 1, 0, 0, 0]])
>

I don't really understand; is this the full Hamiltonian?

> spin_block = tinyarray.array([[0, 0, 0, 1, 0, 0],
>
>      [0, 0, 0, 0, 1, 0],
>
>      [0, 0, 0, 0, 0, 1],
>
>  [-1, 0, 0, 0, 0, 0],
>
>      [0, -1, 0, 0, 0, 0],
>
>      [0, 0, -1, 0, 0, 0]])
>
> And we add “conservation_law= - spin_block” in  “lead=kwant.Builder()”
>
>  
>
> The program can not work, and prompt error:
>
>  
>
> IndexError: index 2 is out of bounds for axis 0 with size 2
>

There is not enough information for me to be able to help you. Please
attach a complete Kwant script.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Wavefunction for 2D NISIN

[Joseph Weston]

Hi,


>  
> I still don’t understand exactly how I can implement this. From what I
> can understand you choose the energy not the eigenstate for the
> density function.


I am unclear what you mean here. In your code you do the following:


    ham_mat = sys.hamiltonian_submatrix(sparse=True, params = params)
    evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=20, sigma=0))

    kwant.plotter.map(sys, np.abs(evecs[:, 9])**2


You take the hamiltonian of the scattering region, diagonalize it, and
then plot *a single eigenvector*. Your problem was that because you have
>1 degree of freedom per site you cannot simply plot the eigenvector
as-is. You need to first compute the density per site for this
eigenvector, which you can do using the operator module:


    rho = kwant.operator.Density(np.kron(-s_z, s_0))  # holes count +1
to the density, and electrons -1

    kwant.plotter.map(sys, rho(evecs[:, 9))


Note that you call the density operator with a single eigenvector.


Independently of the above it's not clear to me that you should be
diagonalizing just the scattering region Hamiltonian, given that your
system has leads. It seems to me that you should be calculating the
scattering states and the density of those.


Happy Kwanting,


Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Difference between hop[0].tag and hop[1].tag

[Joseph Weston]

Hello,

>
> What is the difference between   hop[0].tag   and hop[1].tag
>
hop[0].tag is the tag of the first site in the hopping and hop[1].tag is
the tag of the second site in the hopping.


>
> My second question 
> the onsite fuction for example: def onsite(site, V):return V
>
> why it depends on site where the shape functions for example 
> def circle(pos): rsq = pos[0] ** 2 + pos[1] ** 2
> depends on pos

Because  when creating a shape in realspace you typically only care
about the position, whereas your onsite matrix elements could
potentially depend on other things (e.g. the lattice that the site is from)


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] 3D annulus cylinder

[Joseph Weston]

Hi,


> Dear KWANT developers,
>
> I want to generate a 3D annulus cylinder with hopping in all the three
> directions(periodic in x-direction) like the following figure. Please
> suggest me the way to generate such type of geometry. If possible then
> provide an example code for that for better understanding.
>
> image.png


You can use a 'shape' function as shown in Kwant tutorial 2.3:
https://kwant-project.org/doc/1/tutorial/spin_potential_shape#nontrivial-shapes.
The code will be almost exactly the same as the ring example in the
tutorial, except that you will have 3 coordinates instead of 2.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Questions about conductance

[Joseph Weston]

HI Dominic,

>
> I am using Kwant to stimulate 2D nanowire (square lattice) with two
> barriers (quite similar to NS junction described in Section 2.6 in the
> tutorial but topological superconductor is sandwiched between two
> barriers) and I am measuring the conductance across it using
> "smatrix.transmission(1, 0)". Somehow in a topological phase, I cannot
> get zero bias peak to maintain is something wrong with my code?
> (Checked with spectrum)


You didn't post the code so it's difficult for me to say what the
problem may be. I suspect that you may be calculating the conductance
incorrectly.

You said you used 'smatrix.transmission(1, 0)', however in the presence
of the superconductor this transmission will not give you the
conductance. I think that you should calculate the conductance in the
same way as is shown in tutorial 2.6. Remember to correctly set the
'conservation_law' and 'particle_hole' parameters when constructing the
(non-superconducting) leadsso that you can separate out the electron and
hole subblocks of the scattering matrix.

Post back if something is not clear (or to let us know if this solved
your problem).


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] honycomb lattice with hole

[Joseph Weston]

> sym = kwant.TranslationalSymmetry(latt.vec((-1,0)))
>     
> sym.add_site_family(latt.sublattices[0], other_vectors=[(-1, 2)])
> sym.add_site_family(latt.sublattices[1], other_vectors=[(-1, 2)])
> lead = kwant.Builder(sym)
>  
>  
>

>
> Making a hole on the scatter reign is ok. I did the same procedure for
> lead to make hole on it but at the end there isn't  any hole on the
> leads (it is like a rectangular). That is very kind of you if let me
> know what is the problem? how can I make it?
>

You haven't made the translational symmetry in the leads big enough.
Think about how big the unit cell of the leads is, and how big the hole
you are trying to make is.


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] honycomb lattice with hole

[Joseph Weston]

Hi,


> I’m new in kwant and I want to use this package to calculate transport 
> properties of some graphene-liked materials with multi orbital per atom. In 
> order to define hopping in the form of matrix, we need to determine the 
> relative coordinate of neighbor sites. If we build the system using 
> latt = kwant.lattice.honeycomb
>
> Could anyone help me to show how can I make a hole on the center of
> main region and leadS (I would like to make a graphene nanoribbon
> which have holes)?like in attached fiqure. Please find attached file.
>

Making a hole in a system is shown in the graphene tutorial [1].

Happy Kwanting,

Joe

[1]: https://kwant-project.org/doc/1/tutorial/graphene



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Staggered Sublattice Potential

[Joseph Weston]

Hi,


> Respected Sir,
>
>                      I want to include the staggered sublattice
> potential in Kwant. The problem is, it is an onsite term which takes
> positive and negative values to two different sublattices of Graphene.
> how will I include this onsite potential?
>

I would recommend following the Kwant tutorial. Once you have understood
the main concepts of building Kwant systems, e.g. sites on
(sub-)lattices, you will be able to answer this question yourself.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] A question about the representation of the wavefunction

[Joseph Weston]

Hi Adel,

I read through this thread and I'm still not really sure what is being
proposed. In the first email in the thread it seems as if you are
proposing a wrapper for wavefunctions so that they can be indexed by
site. It seems that you want to be able to write:


   psi(i).dagger() * syst.hamiltonian(i, j, params) * psi(j)


This is a fine proposal. I could imagine that a simple implementation
for such a wrapper would not need to use any exotic storage formats,
C-level accessors or anything: finalized systems already contains an
attribute 'site_ranges' that allows for constant-time access to the
orbital offset of a site. I actually proposed such a wrapper some time
ago, but my proposal was rejected [1] because it did not give vectorized
access to the wavefunction elements. It is not obvious to me how to give
vectorized access without copying.

Then after this proposal the thread seems to start discussing
_LocalOperator (which is an internal detail that was never meant to be
extended or user-facing) and lots of details about implementation
strategies in C++ for *something*, but it is not clear to me what the
end goal is. Perhaps you want to be able to more easily express
arbitrary operators? Or are you proposing a refactoring of the kwant
operators module to make it more readable? A clarification would be very
useful for me.

Thanks,

Joe

[1]: https://gitlab.kwant-project.org/kwant/kwant/merge_requests/47



On 5/14/19 3:14 PM, Adel Kara Slimane wrote:

> Hello Christoph,
>
> I have a few additions to make above my previous mail, I hope I am not
> making you waste your time maybe with trivial and unuseful suggestions.
>
> I have found a C++ library, Armadillo ,
> that can do the trick for implementing a fast accesser, with built-in
> matrix-vector products.
>
> Let's assume we have the wavefunction and the hamiltonian on flat
> contiguous arrays. This library has Row, Column, Matrix and
> SparseMatrix object implementations, with inner functions like
> hermitian conjugate or transpose, and with operations between them
> like sum and multiplication. The good thing about it is that one
> initialise such objects with pointers to existing memory, without
> copying, here are the constructors that one could use:
>
> For vectors: 
>
> |vec(ptr_aux_mem, number_of_elements, copy_aux_mem = true, strict
> = false)| 
>
> For matrices :
>
> |mat(ptr_aux_mem, n_rows, n_cols, copy_aux_mem = true, strict =
> false)| 
>
> So one can wrap such objects in Cython and stack allocate them with
> the correct memory pointers and sizes then use them to write the
> wanted operator expression in a compact way. In terms of speed, I
> expect it to be nearly as fast as a "by-hand" impl! ementation, since
> the expression will be cythonised/compiled, thus in the end with
> similar nested for-loops after compilation.
>
> If I understand well, matrices, for example the hamiltonian, don't
> have really a low level representation since "params" need to be
> provided first to create a flat matrix of only complex values, and
> that would be the point of "_eval_hamiltonian()" in the
> "_LocalOperator" class, it returns a BlockSparseMatrix which is
> basically a flat array containing all the matrices of all the needed
> hoppings, given in the "where" list. I have a question about that:
>
> Would it be as fast to compute the hamiltonian little by little while
> calculating the operator's value on each hopping of the "where" list
> ?  To show what I mean exactly, I have a code example
>  (I
> couldn't test it yet unfortunately) using the "tinyarray" pythong
> package. The idea would then be to use Aramadillo objects instead of
> tinyarray, and have the code compilable. This approach would have the
> benefit of not using "_eval_hamiltonian()" nor the "BlockSparseMatrix"
> class, but will need the extra Armadillo wrapping code, and mapping
> the return of "syst.hamiltonian" to a Matrix object.
>
>
> Thank you,
>
> Adel KARA SLIMANE
>
> On Mon, 13 May, 2019 at 4:38 PM, Adel Kara Slimane
>  wrote:
>> Hello Christoph,
>>
>> Thank you for your! detailed answer.
>>
>> I convinced myself with a small C++ benchmark code that summing over
>> the complex values of a vector>> array is more
>> than twice slower than summing over a flat, same size, dynamically
>> allocated array (I attached the code to this mail). About the memory
>> usage, I don't understand why it would take signifcantly more space,
>> but maybe because I first suggested an array of Python objects
>> instead of compiled objects, which take more space ?
>>
>> I agree and I am interested by coding the solution you are
>> suggesting, I would like to suggest further developpement to it:
>> instead of returning a simple list of 2D arrays without copying
>> anything, we use an "accesser" object 

Re: [Kwant] Error with params SimpleNamespace

[Joseph Weston]

Hi,


> I’ve been running into the same error message over and over when
> trying to pass a SimpleNamespace to the function
> hamiltonian_submatrix(). I am using Python 3.7.3 and Kwant 1.4.1.
> Could it be that this is due to a too recent Python version?
>
>
>  
>
> pars = SimpleNamespace(t=1, mu=-0.1)
>
> hamil = syst.hamiltonian_submatrix(params=pars)
>

When passing parameters using 'params' you need to pass a *dictionary*.
This is in the docstring of 'hamiltonian_submatrix' and all other kwant
functions that accept 'params'. Using SimpleNamespace was a trick that
people sometimes employed before it was possible to pass named
parameters via 'params'.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Peierls phase substitution on Tight-Binding Hamiltonian

[Joseph Weston]

> And also, when I am plotting the spectrum in varing magnetic field
> with the hamiltonian as a sparse matrix. I am getting an value error.
> *ValueError: x and y must have same first dimension, but have shapes
> (100,) and (1, 15) *
> I don't know what does it mean?
>
> *def plot_spectrum(syst, Bfields):
>     energies = []
>     for B in Bfields:
>    
>     ham_mat = syst.hamiltonian_submatrix(params=dict(B=B),
> sparse=True)
>     ev = sla.eigsh(ham_mat.tocsc(), k=15, sigma=0,
>     return_eigenvectors=False)
>     energies.append(ev)
>     pyplot.figure()
>     pyplot.plot(Bfields, energies)
>     pyplot.xlabel("magnetic field [arbitrary units]")
>     pyplot.ylabel("energy [t]")
>     pyplot.show()*
> *plot_spectrum(syst, [iB * 0.002 for iB in range(100)])*


The error message seems pretty clear to me; I am not sure where the
confusion is, especially as the traceback will point to the exact line
in the code where the error occurs. If I were you I would first take the
time to really look at the code that you are writing and understanding
what it is doing. Remember that this is not a general mailing list for
coding help: it is specifically for asking questions about Kwant.


Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Peierls phase substitution on Tight-Binding Hamiltonian

[Joseph Weston]

> Dear Joseph,
>
> As per your instructions I have build the system, but the code is not
> giving the desired result.
> Here is the code please have a look on it.


'y' should be the y position of the sites of the hopping, not an input
parameter to the system construction:

    def hoppingz(site0, site1):

 y = site0.pos[1]

 ...


Happy Kwanting,

Joe





signature.asc
Description: OpenPGP digital signature

Re: [Kwant] STM tip

[Joseph Weston]

> I want to use KPM to calculate LDOS vs. energy for a given position.
> Apparently, I can do it using 'where' in kwant.operator.Density, but
> the input must be a crystal position. I wonder if there is any form to
> find the closest crystal position for a given position vector that
> does not necessarily belongs to the crystal.
Lattices have a method 'closest' that does that.

> Alternatively, can I somehow integrate kwant.operator.Density in a
> specific region and not over the entire system?

Yes, if you provide a 'where' and 'sum=True' then the density will be
summed over the sites specified by 'where'. Note that 'where' may be a
function , so you don't have to specify the sites explicitly.

Happy Kwanting,

Joe





signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Peierls phase substitution on Tight-Binding Hamiltonian

[Joseph Weston]

>
> Then as you said I have to multiply the hopping in z direction by
> exp(2 * pi * 1j *  B * a* y). But how can I build a 3D system from
> here. I know that 2D systems can be written as *lat =
> kwant.lattice.square(a), *but I don't know how should I proceed to
> bulid a lattice in 3D using the above onsite and hoppings? Please help
> me regarding this.


The tutorials explain how to make lattices in arbitrary numbers of
dimensions. How about making a cubic lattice in 3D?


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Peierls phase substitution on Tight-Binding Hamiltonian

[Joseph Weston]

Hi,

On 5/5/19 1:08 PM, Naveen Yadav wrote:
> Dear KWANT developers,
>
> I have an Hamiltonian
> *H(k) = tx*σx* sin kx + ty*σy*sin ky + mk*σz + λ*σ0 *sin kz,*
> *
> *
> *mk = tz(cos β − cos kz) + t'(2 − cos kx − cos ky)*
> **
>
> I am trying to introduce magnetic field in *x-direction* for a gauge
> *A = (0, 0, By),*
> Can I make a *template* for this hamiltonian using KWANT discretizer
> or I have to discretize it by hand to get the onsite and hopping matrix?


I already answered this question in the previous thread. The Hamiltonian
is essentially already discretized; you just have to identify which
terms are onsites and which are hoppings. See the previous thread.


> And then how to add the peierls phase to the hopping matrix?

For the gauge you have chosen the peierls phase will be non-zero only on
hoppings in the z direction. You can specify these hoppings as a
function and multiply the hopping by exp(2 * pi  * 1j * B * a * y) where
"a" is the lattice discretization length and "y" is the y position of
the sites at either end of the hopping (will be the same for both sites,
as we only apply this phase to hoppings in the z-direction).


> Then I want to plot the spectrum as a function of k_z, with Wx =10 and
> Wy=50

The resulting Hamiltonian will be a function of k_z. You can use
kwant.plotter.spectrum to plot the spectrum as a function of k_z.


Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] [KWANT] Diagonalization of hamiltonian

[Joseph Weston]

>
> I have build the system(discretized in x and y direction). The actual
> problem is that, " How to add magnetic field term(as given by the
> auther l_m = 4.5) to this system?" These bands are in presence of
> magnetic field.


Well you'll have to look in the paper to answer that question; you're in
a better position to answer that than me. Maybe a Peierls phase?




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Regarding the NNN term in bilayer graphene across the layer

[Joseph Weston]

>   I have constructed the lattice with third dimension
> Z0 using kwant.lattice.general. As the system is giving a 3d plot, I
> am unable to visualize the hopping between nearest neighbors and so on
> which make me confused whether the hopping kind I had given is right
> or wrong. Is there any way to visualize those hoppings in the 3d plot
> (as we can see in the 2d plot).
>

Ah, that's a good point. Perhaps you could make the 2 lattices in 2D but
put an offset in the x-y positions of the sites in one of the layers of
the bilayer. This way you can still visualize the hoppings using
kwant.plot and the sites occupy distinct positions in space. Because we
are limited by matplotlib's poor 3D rendering we cannot display the
hoppings in 3D plots.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Regarding the NNN term in bilayer graphene across the layer

[Joseph Weston]

Hi,


>  I am trying to put next nearest neighbor hopping terms across the
> layer in bilayer graphene. In that case, how can I get those hopping
> kind directions for the NNN hopping for the interlayer or can I use
> sys [ lat.neighbor(3)] or sys[ lat.neighbor(4)]?
>
The answer will surely depend on how you have constructed your bilayer
graphene. Kwant does not have a built-in way of constructing a specific
bilayer honeycomb lattice, so you will have needed to construct this
yourself from 2 honeycomb lattices. If you don't mind the sites from the
2 lattices living at the same realspace position then the following
could work:


    lower_layer = kwant.lattice.honeycomb(name='lower-layer')

    upper_layer = kwant.lattice.honeycomb(name='upper-layer')

    # Loop over the A and B sublattices from the upper and lower layers

    for upper, lower in zip(upper_layer.sublattices,
lower_layer.sublattices):

        inter_layer_hoppings = kwant.builder.HoppingKind((0, 0), upper,
lower)

        syst[inter_layer_hoppings] = ...


If you want to introduce a third dimension, so that the lower and upper
layer sites are actually offset in realspace then you can do exactly the
same thing, except that you'll have to construct the honeycomb lattices
using 'kwant.lattice.general', giving the primitive vectors in 3D (zero
z component) and the basis in 3D also.


Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] [KWANT] Diagonalization of hamiltonian

[Joseph Weston]

Hi,


> I want to diagonalize the model hamiltonian containing sine and cosine
> functions with momentum operators as their argument.
>
> H(k) = tx*σx* sin kx + ty*σy*sin ky + mk*σz + λ*σ0 *sin kz,
>
> mk = tz(cos β − cos kz) + t'(2 − cos kx − cos ky)

If you just want to diagonalize a 2x2 H(k) then you don't even need
Kwant; you can just make a function that gives you H(k) given a k
vector, and then diagonalize the result using scipy.linalg.eigs.


> Then I want to plot the energy dispersion as a function of kz in the
> presence of perpendicular magnetic field. Here perpendicular direction
> is x.

Does the model already contain the terms for a magnetic field in the x
direction? The answer will depend on what effects you are taking into
account: do you want the orbital component, or just the action on the
spin degree of freedom? The answer will depend on what you are trying to
model, and this is a question that you will need to answer for yourself.


Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] discretised Hamiltonian

[Joseph Weston]

Hi,


> Dear Developer,
>          The discretised Hamiltonan shown in Kwant release version
> 1.4.0 in Article# 2.2.1 is not being understood when I try to attach
> significance of "ket" and "bra" state to the symbols e.g.   │i+1,j>
>  and |i+1,j>. Such operators can be used as a basis for
writing down the Hamiltonian as matrix.


>  The discretised  1D -Hamiltonian is given like :
> H= [-tΨ_n-1 -2tΨn +tΨ_n+1 ] /a^2  
> "t" being hopping integral and "a' the lattice constant.


To me this looks more like you're applying the Hamiltonian to a
wavefunction.


> In the article, │i,j> is the positional state. Then what signifies the
> "bra"? I mean,  The question would appear very simple to understand that I know. Yet,
> I am confused. 
> Please, help me to rightly think of the meaning of the Hamiltonian
> expression.
> Happy Kwanting and looking forward to fix my correct understanding.


The bra is a linear functional on the space of states. Specifically,
. This sounds like a misunderstanding of some underlying
concepts; I would recommend reading a text on the foundations of quantum
mechanics.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

[Kwant] Kwant tutorial notebook

[Joseph Weston]

Dear all,


I recently gave a tutorial at the University of Maryland about quantum
transport and more specifically using Kwant in that context. The
tutorial also contained a relatively in-depth discussion of the linear
algebra that Kwant performs internally.

Because this is of general interest to the community I have made  the
slides public on GitHub [1] , and also in executable form on Binder [2].
The slides are in the form of a Jupyter notebook, and released under a
Creative Commons license, so you can download run and modify to your
heart's content!


Happy Kwanting,

Joe


[1]: https://github.com/jbweston/maryland-kwant-tutorial/

[2]:
https://mybinder.org/v2/gh/jbweston/maryland-kwant-tutorial/master?filepath=index.ipynb




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] kwant installation via conda in windows 10

[Joseph Weston]

Hi Xiaoyan,

> In the newest release of kwant (Ver 1.4.0), it has been announced that
> windows installation can be done via Conda. However, there is no guide
> about this in the “installation section” page. It seems that
> “installation section” is old and not matching the new release. Any
> step-by-step guide on conda installation in windows? Thanks.
>
I just created a merge request for the website for this; thanks for
reminding us!

https://gitlab.kwant-project.org/kwant/website/merge_requests/29


Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Attaching single leads - quantum transport

[Joseph Weston]

Hi Nuwan,


> I am new to Kwant and I am trying to attach a single infinite lead
> from both sides for the case of quantum transmission through a single
> impurity site to get the maximum transmission probability. 
>
> Also, how can we develop the tutorial example of "Transport through a
> Quantum Wire", when attaching 'single infinite leads' to it?

In Kwant you need to treat the "single infinite lead [attached] from
both sides" as 2 separate leads, so just follow the example from the
Kwant tutorial.


Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Transmission matrix and momenta

[Joseph Weston]

Hi Luca


> Dear all,
>
> I am trying to obtain the coefficients of
> the transmission matrix (fine, so far.)
> and to related them to the modes of the leads.

This is covered in the FAQ [1]. You can use 'lead.modes()' to get the
mode wavefunctions, and the scattering matrix entries are ordered in the
same way as the modes. i.e. the columns are ordered as the incoming
modes of the associated leads, and the rows are ordered as the outgoing
modes of the associated leads.


> In particular, I would like to know the momenta of these
> modes (supposing periodic conditions along the directions orthogonal
> to the leads).
> Is it possible ?

I'm not sure why the periodic boundary conditions are relevant?
'lead.modes()' will give you the momentum along the translational
symmetry direction.


Happy Kwanting,

Joe


[1]:
https://kwant-project.org/doc/1/tutorial/faq#how-does-kwant-order-the-propagating-modes-of-a-lead




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Updating Kwant

[Joseph Weston]

Hi Marc,

> I want to implement the latest development of Kwant since I want to
> try to calculate Kubo conductivities with the KPM method (as it
> appears in the development version of the documentation). I tried to
> use Conda to update the packages (conda install -c kwant kwant) but I
> still get the message"AttributeError: module 'kwant.kpm' has no
> attribute 'conductivity'".
>

Unfortunately the development conda builds have been failing for some
time and we didn't get round to fixing it (but Kwant itself still
works). As Christoph says, we're in the process of making a release. You
can track how we're doing in this issue:
https://gitlab.kwant-project.org/kwant/kwant/issues/275. We'll probably
be able to fix the development builds as part of/just after the release,
but by then you'll be able to get a copy of Kwant 1.4 anyway.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Berry phase of a system periodic in 1 direction

[Joseph Weston]

Hi,

On 2/15/19 6:42 PM, Srilok Srinivasan wrote:
> Dear Kwant users and developers,
>
> Is there a way to calculate the integral(< u_nk | d_k(u_nk)>dk) 
> across the 1D BZ of a system periodic in one direction. Here u_nk is
> the periodic part of the Bloch wave function corresponding to band n.
> I am trying to compute this for a simple system of graphene nanoribbon.
>
> I am trying compute the gauge independent value of the above integral
> using a "1D Wilsoon loop" where the integral is approximately equal to
> the product of   for all k across the BZ.  I defined
> the unit cell of the nanoribbon as a lead and tried to get the u_nk as
> the propagating modes using leads.modes() method (I have pasted the
> code below) but the resolution in k is too large - only 4 momenta
> across the BZ.

'lead.modes' computes all the modes at a given energy (0 by default).
This is probably not what you want. What you really want is to compute
the eigenvectors of H(k) for different k. 'kwant.physics.Bands' does
this internally, but in Kwant 1.3 it doesn't give you access to the
eigenvectors (just an oversight on our part).

We're in the process of releasing Kwant 1.4, where Bands has been
modified to also return the eigenvectors, so you could wait to use that
(you can track our progress here:
https://gitlab.kwant-project.org/kwant/kwant/issues/275). Alternatively
you can just construct H(k) = h + v exp(ik) + v^dagger exp(-ik) where h
is 'lead.cell_hamiltonian()' and v is 'lead.inter_cell_hopping()' and
diagonalize it yourself for different values of k.

Having said that, you need to be careful. To the best of my knowledge
numerical eigensolvers will typically give you the eigenvectors in
different gauges at different k values (i.e. there will be an arbitrary
phase factor), so  naively calculating the derivative using some finite
difference scheme might give nonsense (I'm not sure about this, though).


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Problem with semimetals

[Joseph Weston]

Hi Luca,



> I have tried to implement a couple of examples of semimetals with Kwant,
> where the lattices (scattering regions) are taken square or cubic,
> and the leads attached on entire faces of the lattices.
> I assumed both open and periodic (in the directions orthogonal to the leads)
> boundaries conditions for the lattices and for the leads.

I'm not sure what exactly you mean here. Could you post a code snippet
that illustrates what precisely you mean?

> Discussing at a conference,
> I got aware vaguely about possible problems
> that may arise when one wants to
> discretize in real space an Hamiltonian originally written in
> the space of momenta (also continuous).
> This is exaclty what I have done always;
> notice that I performed the discretization by hand, without
> exploiting the routine in Kwant.
>
> Do you know about similar issues with semimetals ?
> Do you have some advices ?

So what you're saying is that you start with a continuum Hamiltonian and
discretize it, correct?

For certain classes of models this can certainly lead to unphysical
states, even at energies where the continuum model is valid.

For example, if we take a continuum model for graphene near a single
dirac point so H∝|k|, discretize it with finite differences and
diagonalize the resulting discrete Hamiltonian, we will see that the
spectrum is proportional to sin(k)! This means that even at low energies
(where the continuum model is valid) there are states close to k=π,
which are *not* present in the continuum model, and that were introduced
entirely as a result of the discretization.

The general term for this is the "Fermion doubling problem". You should
be able to tell if this is your issue by comparing the spectrum of the
Kwant system with the spectrum of the continuum model from which it is
derived. If you see that the "bands bend down" to give spurious states
even at low energies then this could indicate that you have a Fermion
doubling problem.

You have to be more careful with the discretization to avoid the Fermion
doubling problem, for example section B of
https://arxiv.org/abs/0810.4787. I am not an expert in this field, so if
it turns out that your problem *is* Fermion doubling, I would look into
the relevant literature for that.


Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Datta’s ballistic transport formalism vs KWANT

[Joseph Weston]

Hi Ran,

> I tried to compare KWANT’s results for transmission with Datta’s
> ballistic transport formalism where total transmission is written as
>
> Ttot=T(E)M(E)
>
> Here Datta takes T(E)=1 for ballistic transport (please see: J. Appl.
> Phys. 105, 034506, 2009) and M(E) is the number of modes in transverse
> direction. When I compared KWANT's results with Datta’s expression,
> for the system given in “quantum_wire_revisited.py”, I found different
> results (please see the attached figure where I tried to put every
> relevant thing in the calculation). Since the reflectance is zero for
> that system and so transmission is 1 for each mode, shouldn’t it give
> the same results with Datta’s transmission expression?
>

Nice question!

Looking at your results it seems that the energies at which new modes
open is shifted with respect to Datta's result.

I believe that this is simply due to the fact that your discretization
is not fine enough. Datta's result is valid in the continuum limit,
whereas the Kwant simulation (in the case presented) uses a
finite-difference discretization to render the problem discrete. If you
decrease the 'a' parameter, you should see the discrepancy between the
two result decrease.


Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Install Kwant on Win 10

[Joseph Weston]

Hi,



>
> Thanks for reply. Actually I installed python 3.6 and can do python
> programming on my PC.
>
> image.png
>
> Do you have any idea about this?
>

Probably Python is not in your PATH. Did you exactly follow the
instructions on the Kwant website?

It seems that step 4
here:https://kwant-project.org/install#microsoft-windows should add Python

to your PATH so that you can type 'python'. rather than
'C:\Python36\python.exe' every time.


In any case this problem is with your Python installation, and is not
specific to Kwant.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] complex conjugate operator

[Joseph Weston]

Hi 张金龙,

Thanks for posting to the mailing list.

>  Recently, I want to calculate the conductance of Superconductor
> junction.
>
> Let me take the p+ip topological superconductor as the example. Its
> Hamiltonian is like
> with \sigma, \tau is the Pauli matrix for spin and particle-hole.
>
> It have particle-hole symmetry: 
> to g
> uarantee
>  
> .
> **
> My question is how to write down the K (complex conjugate) operator?


I am not sure what you mean by "write down". You cannot write it down as
a matrix; it's defined by its action on wavefunctions.

I'm not even sure that I understand your question. What are you actually
trying to do? Perhaps reading the tutorial on superconductivity will
help: https://kwant-project.org/doc/1/tutorial/superconductors. The
example in that tutorial shows how to tell a Builder that your system
has particle-hole symmetry.


Happy Kwanting,

Joe

Re: [Kwant] Drawing edge states in zGNR (using tutorial 2.8)

[Joseph Weston]

Hi Elmo,


> Is there a simple example out there showing how I can build a zigzag
> nanoribbon, just using nearest neighbour hopping, and extract the
> wavefunction and draw said nanoribbon exactly like the figure at the
> end of the 2D part of that tutorial.


Can't you just take the example from the tutorial and then make a
nanoribbon shape instead of a circle shape?


Happy Kwanting,

Joe

Re: [Kwant] semicon

[Joseph Weston]

Thanks for the information!

> Thank you for your quick response. I am using kwant 1.3.3 version and
> semicon 0.1.0 version. I think the error is related to calling
> parameters function which are stored as a tuple. By the way, this is
> the error I encountered:
>
> TypeError Traceback (most recent call last)
> /home/ali/anaconda3/lib/python3.6/site-packages/kwant/builder.py in 
> hamiltonian(self, i, j, params, *args)
> 1995 try:
> -> 1996value = value(self.sites[i], **params)
> 1997 except Exception as exc:
>
>  in onsite(site, Delta_0, E_0, E_v, P, gamma_0, gamma_1, gamma_2, 
> hbar,
> k_y, m_0)
>
> TypeError: __call__() takes 2 positional arguments but 3 were given
>
> The above exception was the direct cause of the following exception:
>
> UserCodeError Traceback (most recent call last)
>  in ()
> 59 # compute the scattering matrix at a given energy
> 60 p = {'k_y': 0, **two_deg_params}
> ---> 61smatrix = kwant.smatrix(syst, energy, params=p)
> 62 
> 63 # compute the transmission probability from lead 0 to
>
> /home/ali/anaconda3/lib/python3.6/site-packages/kwant/solvers/common.py in 
> smatrix(self, sys, energy, args, out_leads, in_leads,
> check_hermiticity, params)
> 370 linsys, lead_info = self._make_linear_sys(syst, in_leads, energy, 
> args,
> 371   check_hermiticity, 
> False,
> --> 372params=params) 373 
> 374 kept_vars = np.concatenate([coords for i, coords in
>
> /home/ali/anaconda3/lib/python3.6/site-packages/kwant/solvers/common.py in 
> _make_linear_sys(self, sys, in_leads, energy, args, check_hermiticity,
> realspace, params)
> 162 lhs, norb = syst.hamiltonian_submatrix(args, sparse=True,
> 163return_norb=True,
> --> 164params=params)[:2] 165 lhs = getattr(lhs, 'to' + 
> self.lhsformat)()
> 166 lhs = lhs - energy * sp.identity(lhs.shape[0], 
> format=self.lhsformat)
>
> kwant/_system.pyx in kwant._system.hamiltonian_submatrix()
>
> /home/ali/anaconda3/lib/python3.6/site-packages/kwant/builder.py in 
> hamiltonian(self, i, j, params, *args)
> 1996 value = value(self.sites[i], **params)
> 1997 except Exception as exc:
> -> 1998_raise_user_error(exc, value)
> 1999 else:
> 2000 try:
>
> /home/ali/anaconda3/lib/python3.6/site-packages/kwant/builder.py in 
> _raise_user_error(exc, func)
> 1883 msg = ('Error occurred in user-supplied value function "{0}".\n'
> 1884'See the upper part of the above backtrace for more 
> information.')
> -> 1885raise UserCodeError(msg.format(func.__name__)) from exc
> 1886 
> 1887 
>
> UserCodeError: Error occurred in user-supplied value function "onsite".
> See the upper part of the above backtrace for more information.


this certainly looks like a bug in semicon. I've managed to reproduce
the bug, but I don't understand why it occurs. I inspected what was
happening using the Python debugger and it seems to me that everything
should work (but clearly it doesn't!).


Unfortunately the author of semicon is no longer working in academia,
and the semicon project is currently without a maintainer (and it's not
even v1.0!), so I am not sure what the best course of action is going
forward.


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] semicon

[Joseph Weston]

Hi Ali,


> I want to find conductance of a heterostructure using semicon library.
> When I just use coords='z', everything is OK. But, when I want to
> change it to two dimensions like 'xz', I encounter an error. Would you
> please tell me how I can fix it? I hope it is not just a
> simple mistake which I am not aware. Your time and consideration are
> greatly appreciated in advance. Please find attached the code
> correspondent to what I mentioned above.


Thanks for posting the script, but you'll need to give more information
before anyone else can help you:


+ what version of Kwant and semicon are you using

+ what error was produced? Please provide a full traceback


Happy Kwanting,


Joe

[Kwant] Change conductance step units

[Joseph Weston]

The message below was originall from Sergio Castillo. For some reason
neither I nor several other Kwant developers received this message from
the mailing list, so I'm reposting it here.


> Hello everyone, I would appreciate any suggestion you could give me to
> solve this:
>
> Im trying to plot the conductance for a two dimensional material,
> however, I want the y-axis to show conductance in (2e^2/h) units
> instead of (e^2/h) units that are shown by default as in the examples.
> It could be something obvious but I can't figure it out how to plot
> the data in that units
>
> Thanks in advance.

Hi Sergio,

You can divide your conductance  by 2, that will change the units from
(e^2/h) to (2e^2/h)

Happy Kwanting,

Joe


signature.asc
Description: OpenPGP digital signature

Re: [Kwant] How to plot the lead band of non hermitian Hamiltonian with Kwant

[Joseph Weston]

Hi,


> In the kwant, I want to plot the lead band of non hermitian
> Hamiltonian and I just add some imaginary terms to the onsite
> term.However, the error "The cell Hamiltonian is not Hermitian" will
> appear. How can I solve the problem? 
>              I am looking forward for your reply and I would
> appreciate it if you could contact me.


Indeed it's not possible to use the built-in band plotter with
non-hermitian Hamiltonians. If your Hamiltonian is non-Hermitian it will
in general have complex eigenvalues, and there's no way to plot this on
a simple line plot.


What exactly are you trying to accomplish? Are you sure that using
non-hermitian leads is the way to achieve this?


Happy Kwanting,


Joe



signature.asc
Description: OpenPGP digital signature

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

[Joseph Weston]

Hi everyone,


> This is very strange; when I download the notebook and run it I indeed
> see the error that you reported, however when I try to run the
> following minimal example it succeeds:
>
>
>> import numpy as np
>> import matplotlib.pyplot
>> import kwant
>>
>> graphene = kwant.lattice.general([[1, 0], [1/2, np.sqrt(3)/2]],  #
>> lattice vectors
>>  [[0, 0], [0, 1/np.sqrt(3)]])  #
>> Coordinates of the sites
>>
>> zigzag_ribbon = kwant.Builder(kwant.TranslationalSymmetry([1, 0]))
>> zigzag_ribbon[graphene.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = 0
>> zigzag_ribbon[graphene.neighbors(1)] = 1
>>
>> print(kwant.__version__)
>> kwant.plot(zigzag_ribbon);
>
> I'll investigate a bit more and report my results.
>

Found it! In Kwant 1.4.0a1 we tweaked Kwant's plotting so that sites
from different families (sublattices) would use the colors from
Matplotlib's color cycle by default. In the March Meeting tutorial
notebooks there is a line:

> %run matplotlib_setup.ipy
that loads some custom matplotlib setup to make the plots look a bit
better in the notebooks. Looking into 'matplotlib_setup.ipy' we see the
following suspicious line:

> mpl.rc('axes', prop_cycle=mpl.cycler('color', ['black']))
which was included to stop all the bandstructure plots being
multicolored (because it's distracting), however because there are fewer
colors than there are site families (2 site families in graphene) the
color cycling in Kwant's plotter broke.


This is a bug in Kwant, and I've opened an issue [1] on our bugtracker.
I'll fix this today and it'll go into the v1.4.0 release (that we will
make sometime between now and christmas).


In the short term you can fix this in a couple of ways:

+ either downgrade to Kwant 1.3.3

+ or change the offending line in 'matplotlib_setup.ipy' to:

    mpl.rc('axes', prop_cycle=mpl.cycler('color', ['black', 'black']))

+ or wait until I push a change to the tutorial notebooks later today
and then re-download them


Thanks for the report, and happy kwanting!


Joe


[1]: https://gitlab.kwant-project.org/kwant/kwant/issues/257




signature.asc
Description: OpenPGP digital signature

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

[Joseph Weston]

Just to respond to Anton:

> 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

They're not plotting a finalized system:

>>> 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());

so this should totally still work (and indeed it does on Kwant 1.3.3, or
when taking a minimal example as I illustrated in my email.


Joe




signature.asc
Description: OpenPGP digital signature

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

[Joseph Weston]

Hi Srilok,


>
> 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
>
> 
> ---
> KeyError  Traceback (most recent call 
> last)
>  in 
> 5 zigzag_ribbon[graphene.neighbors(1)] = 1
> 6 zigzag_ribbon.finalized()
> > 7kwant.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)
> --> 964site_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):
> --> 909spec = [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):
> --> 909spec = [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):
> --> 958return 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)
>

Thanks for the bug report!

This is very strange; when I download the notebook and run it I indeed
see the error that you reported, however when I try to run the following
minimal example it succeeds:


> import numpy as np
> import matplotlib.pyplot
> import kwant
>
> graphene = kwant.lattice.general([[1, 0], [1/2, np.sqrt(3)/2]],  #
> lattice vectors
>  [[0, 0], [0, 1/np.sqrt(3)]])  #
> Coordinates of the sites
>
> zigzag_ribbon = kwant.Builder(kwant.TranslationalSymmetry([1, 0]))
> zigzag_ribbon[graphene.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = 0
> zigzag_ribbon[graphene.neighbors(1)] = 1
>
> print(kwant.__version__)
> kwant.plot(zigzag_ribbon);

I'll investigate a bit more and report my results.


Happy Kwanting,


Joe

Re: [Kwant] Fixing min and max values of colormap using kwant.plot

[Joseph Weston]

Hi,


> Thank you for the quick response. When I try /kwant.plotter.map(sys,
> value, vmin=0, vmax=2) /I get the following error:
> ValueError: Only 2D systems can be plotted this way.
> My system is a 3D system. I would be glad if you could help me out in
> this regard.
>

As the error message says, 'kwant.plotter.map' doesn't work for 3D
systems. Matplotlib does not have very good support for 3D plots in any
case, and it's not clear that there's a general "best way" of plotting a
3D scalar field.

Probably your best bet will be to reduce the dimensionality of what
you're trying to plot.

Failing that, you could try using 'kwant.plotter.interpolate_density'
(available in the current development version of Kwant and in the
soon-to-be-released Kwant 1.4) with the 'ipyvolume' [1] package;
something like:


    import ipyvolume

    ...

    field, box = kwant.plotter.interpolate_density(syst, some_density)

    field = field.squeeze()  # ipyvolume requires a rank-3 array

    ipyvolume.quickvolshow(field)


Note that the above code is incomplete and untested; it's just to show
how you could glue the relevant pieces together.


Happy Kwanting,

Joe



[1]: https://ipyvolume.readthedocs.io/en/latest/

Re: [Kwant] The right way to build up a system?

[Joseph Weston]

Hi Guangze,


> I am trying to plot the density and current density of a chiral edge
> state in a 2D quantum Hall system. I tried to build up the system in
> two ways:
> 1. define two sublattices,
> 2. define one lattice, use matrix form for the on-site potential and
> hoppings

> I expected that these two ways should yield the same result. However,
> the result for the chiral edge state is not the same. In terms of
> density, method 1 gives the correct picture whereas method 2 doesn't.


This is because in the second case you are only plotting half the
degrees of freedom. There are N sites in your system, but the
wavefunction contains 2N elements (2 degrees of freedom per site).

For this reason kwant provides operators that allow you to extract
per-site quantities. For example, to get the charge density you just need


    rho = kwant.operator.Density(syst)

    kwant.plot(syst, site_size=rho(wf), ...)


And to get the spin density you'd need


    rho_s = kwant.operator.Density(syst, onsite=sigma_z)

    kwant.plot(syst, site_size=rho_s(wf), ...)


> In terms of current density, the resulting picture is the same,
> whereas len(current) differs by a factor of 2.

This is because the current operator calculates current between sites.
In your first system there are 2N sites, and in the second case there
are only N sites.





I would finish by saying that in modern Kwant there is not much reason
to split degrees of freedom onto separate lattices. It is almost always
clearer to write the Hamiltonian using matrix-valued onsites and
hoppings, and using the functionality in 'kwant.operator' means that you
rarely need to manually access wavefunction elements.


Happy Kwanting,


Joe

Re: [Kwant] Extract wavefunction for one lattice

[Joseph Weston]

Hi Yuhao


> In the system I use two lattices to represent spin up and down.
>
> lat_u = kwant.lattice.honeycomb(a=1, name='up') 
> lat_d = kwant.lattice.honeycomb(a=1, name='down')
>
> When I obtain the scattering wave functions from
> wf=kwant.wave_function(sys, en), 
> it contains the value for both lat_u and lat_d.


On 11/01/2018 01:33 PM, Abbout Adel wrote:
> Dear Yuhao,
>
> The order of the wavefunction elements is the same as that of the
> sites. Since you know the order of sites from sys.sites the task
> becomes very simple.
>
> In the following example, you can check the family of the sites in
> sys.sites and notice that you have the electrons first and then the
> holes. The wavefunction order becomes straightforward.
>

Adel's answer is correct, however I wanted to point out that there's
little need to separate internal degrees of freedom onto separate
lattices in modern Kwant.

You can define your system much more concisely using matrix-valued
onsites, not to mention you can use more advanced functionality such as
the kwant.operator module to get spatially resolved spin currents etc.

Historically the only reason to separate internal degrees of freedom
(e.g. spin) onto separate lattices was to force Kwant to choose a
certain mode basis  in the leads, e.g. modes with spin up and spin down,
which would mean that the scattering matrix blocks could be more easily
interpreted (e.g. reflection from spin up channel into spin down
channel). In modern Kwant you can declare that your lead satisfies a
certain conservation law (e.g. spin conservation), as illustrated in the
documentation [1], which forces Kwant to use a mode basis that is
compatible with the conservation law.

For example, if your on-site degrees of freedom are (spin up, spin down)
and your lead Hamiltonian conserves spin, then you could declare your
lead like so:

    sigma_z = np.array([[1, 0], [0, -1]])
    lead = kwant.Builder(lead_symmetry, conservation_law=-sigma_z)
    ...

This will guarantee that each of the lead modes has a well-defined spin,
i.e. that a mode wavefunctions are non-zero only on the spin up or spin
down degrees of freedom.

Happy Kwanting,

Joe


[1]: https://kwant-project.org/doc/1/tutorial/superconductors

Re: [Kwant] Units of density

[Joseph Weston]

Hi,


> I've found in other threads in the mailing list that the units of
> current is for example (unit of charge)/(hbar/unit of energy)
> (https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01100.html).
> Also, the local density of states has units of energy/volume
> (https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00169.html​).
>
>

The units of local density of states is rather "per energy per volume"
(this is what is written in the linked message) not "energy per volume".
Though

> My question is, what is the units of the output of the density
> operator? Is it energy/volume as well?
>

It is "per energy per volume", as is the local density of states.

If you sum the output of a kwant.operator.Density for all scattering
states at a given energy and divide the result by 2pi, it will be
identical (up to numerical precision) to the output of kwant.ldos. I've
attached a script that illustrates this.

> I ask this because intuitively I view it as the square of the
> wavefunction, but it gives me values larger than 1 for each site when
> there is only 1 mode involved (see attached picture) ​so it is not
> just the probability of findinge the electron at that site because
> this should be maximum 1. I have also noticed that the values I get in
> the colorbar depend on the value of my hopping (e.g. case of
> graphene), but overall I'm not so sure of the units.
>

The scattering wavefunctions are not normalized over the scattering
region, so if you sum the absolute square of the wavefunction you will
not obtain 1.  The lead modes are normalized such that they carry unit
current, and the scattering wavefunctions are thus normalized in a way
that is commensurate with this normalization of the lead modes.

In the attached script I also show that the norm of the scattering
wavefunction over the scattering region is not 1.

Happy Kwanting,

Joe
import string
import numpy as np
import kwant

def onsite(site, salt):
return kwant.digest.uniform(site.tag, salt=salt)


def hopping(a, b, salt):
c = bytes(a.tag) + bytes(b.tag)
return complex(*kwant.digest.uniform2(c, salt=salt))


lat = kwant.lattice.square(norbs=1)

syst = kwant.Builder()
syst[(lat(i, j) for i in range(10) for j in range(5))] = onsite
syst[lat.neighbors()] = hopping

lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
lead[(lat(0, j) for j in range(5))] = 4
lead[lat.neighbors()] = -1

syst.attach_lead(lead)
syst.attach_lead(lead.reversed())

syst = syst.finalized()


# Loop over a bunch of realizations of the randomness

params = dict(salt='')
energy = 0.5

for params['salt'] in string.ascii_lowercase:
wf = kwant.wave_function(syst, energy=energy, params=params)
psis = np.array([psi for i in (0, 1) for psi in wf(i)])
norms = np.sum(np.square(np.abs(psis)), axis=1)  # sum over all orbitals
# some norm may be 1 by coincidence, but this is unlikely.
assert not any(np.isclose(norms, 1))

rho = kwant.operator.Density(syst).bind(params=params)
op = sum(rho(psi) for psi in psis) / (2 * np.pi)

ldos = kwant.ldos(syst, energy=energy, params=params)

np.testing.assert_array_almost_equal(op, ldos)


print('all equal!')

Re: [Kwant] issues with Current calculation

[Joseph Weston]

Hi Sergey,

> Also, I found that current asymmetry was indeed there and it was
> caused by numerical errors: the precision of the current calculation
> is higher near the lead for which it is calculated, while the error
> may grow significant near other lead(s).  In my case,  at some
> energies it is symmetric,  at some - not.

This is very surprising; do you think that this could be a bug in Kwant?
Can you share a simple script that reproduces this issue (the simpler
the better; that way it's easy to debug)?

Thanks,

Joe

Re: [Kwant] issues with Current calculation

[Joseph Weston]

Hi Sergey,


> 1)  This code gives a strange error:
>
> The debugged program raised the exception unhandled TypeError
> "streamplot() got an unexpected keyword argument 'start_points'"
> File:
> /home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py,
> Line: 2164
>
> which is strange since start_points  is in the streamplot documentation

Ah, in the documentation [1] it does indeed mention "streamplot", but we
are actually referring to 'kwant.plotter.streamplot', not matplotlib's
one. I will modify the documentation so that the distinction is clearer.
Our version of streamplot includes things like including the background
color etc.

We don't currently provide a way to customize all the options of the
underlying calls to matplotlib, because this would be excessively
complicated. Luckily it is pretty simple to write your own wrapper that
you can customize:

    def custom_current(syst, current):
    field, box = kwant.plotter.interpolate_current(syst, current,
...)  # your custom args here

    # make grid for streamplot
    X = np.linspace(*box[0], num=field.shape[0])
    Y = np.linspace(*box[1], num=field.shape[1])

    pyplot.streamplot(X, Y, field[:, :, 0], field[:, :, 1], ...)  #
your custom args here


Make sure you post back to the mailing list if you find some heuristic
that makes decent plots!



> 2)  How to make several plots be generated without manual closure of
> the previous one?  My command     pyplot.show(block=False)  does not
> seem to work in this case.

You can pass a matplotlib axis to 'current' via the 'ax' parameter. This
will make the plotter draw into the axes that you provide.

You'll need to search on stackoverflow or look in the matplotlib docs to
find out how to do exactly what you want, because even though you can
totally define several plots, I'm not sure matplotlib can display them
simultaneously.

Happy Kwanting,

Joe


[1]:
https://kwant-project.org/doc/1/reference/generated/kwant.plotter.current

Re: [Kwant] issues with Current calculation

[Joseph Weston]

Good morning,


> Hi Joe,  Thank you, now it works!

Nice!

>   But I am a bit confused about the physical meaning of the current,
> created by modes in lead n.   How to plot the real current between the
> two leads?  Imagine,  as the simplest case, that my geometry is
> inversion-symmetric,  so that lead 0 goes into lead 1 under
> inversion.  Then I expect the current to be inversion-symmetric as well.

Summing the current contributions for *all* the scattering states
originating from lead L at an energy E gives you the current you would
measure if the fermi levels in all the leads are at energy E and you
apply an infinitesimal voltage to lead L.

For example if we want to see the current profile if we add an
infinitesimal voltage dV to lead 0:

    wfs = kwant.wave_function(syst, energy=E, params=...)
    J_0 = dV * sum(J(psi) for psi in wfs(0))

This is analogous to the transmission obtained by kwant.smatrix:

    smatrix = kwant.smatrix(syst, energy=E, params=...)
    I_10 = dV * smatrix.transmission(1, 0)

I_10 is the current we would measure in lead 1 after applying dV to lead 0.

In the above I have elided the e^2/h factor for brevity, but hopefully I
have been clear enough to get the point across.

Hope that helps,

Joe

Re: [Kwant] issues with Current calculation

[Joseph Weston]

You need to provide the parameters to the current operator when you call
it, because it needs to evaluate the Hamiltonian to calculate the
current (IIRC this is in the docs).


Happy Kwanting,


Joe


On 10/18/2018 06:16 PM, Sergey Slizovskiy wrote:
> Thank you, Joseph,
>   It could be nice to have this in tutorials.
> In my case,  I still get a error
>  Exception "unhandled kwant._common.UserCodeError"
> Error occurred in user-supplied value function "hopping".
>
> Although, my hoppings are pretty innocent  and conductance calculation
> works fine
>
> On 18/10/18 17:04, Joseph Weston wrote:
>> currents = [J_0(p) for p in psi]
>
> Thanks,
>
> Sergey
>

Re: [Kwant] issues with Current calculation

[Joseph Weston]

Hi,

>   I am trying new Kwant 1.3 features on a very simple example of
> graphene monolayer, and, apart from a standard well-working code I add:
>
> for x in sys.sites(): x.family.norbs=1 for x in leadxp.sites():
> x.family.norbs=1 for x in leadxm.sites(): x.family.norbs=1 (to avoid
> norbs not defined error)
>
> ...
>
> J_0 = kwant.operator.Current(sys) wf =
> kwant.solvers.default.wave_function(sys, energy=0.01,
> args=[phi,chempot,0]) psi=wf(1) current = J_0(psi)  
>

'psi' is a bunch of wavefunctions, but the operators only work on one
wavefunction at a time, so you'd have to :

    currents = [J_0(p) for p in psi]

to get the current for each wavefunction.

This is actually a pretty common thing to want to do, so it would make
sense for the operators to be able to work in a "vectorized" manner.

I'll open an enhancement issue on the bugtracker.

Happy Kwanting,

Joe

Re: [Kwant] Attaching 3D Leads to a 3D System

[Joseph Weston]

Hi,


>
> I am trying to attach a 3D lead to a 3D system that I have built, but
> am facing troubles in the same. I am aware of another thread
> (https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2014-May/000125.html)
> highlighting my question, but having gone through that thread I was
> still unable to successfully attach a lead to my system.
>
> Ideally, I want to attach gold leads to my system but I started by
> attaching a simple cubic lattice (by breaking it into 2 different
> lattices). Following is my code.
>

You define your scattering region over what you call 'lat0' and 'lat1',
which are 2D lattices *embedded in 3D space*.

You then create a new lattice 'lat0_lead' using 'kwant.lattice.square',
which creates a 2D lattice *embedded in 2D space*.

You then add sites from 'lat0_lead'  to the scattering region, so the
scattering region now contains sites that are embedded in 2D realspace
and 3D realspace. Kwant is totally fine with this, but the plotter does
not know how to plot this, so an exception is raised.

You can easily get around this by defining 'lat0_lead' to be a 2D
lattice embedded in 3D:

    # 2 orthogonal lattice vectors in 3D space
    lat0_lead = kwant.lattice.general([(dy, 0, 0), (0, dy, 0)])
    sym_lead0 = kwant.TranslationalSymmetry((dy, 0, 0))


The exception raised is:

> ValueError: Input has irregular shape.

This is a bit cryptic, and I've opened an issue [1] against this
usability bug,



Happy Kwanting,

Joe

P.S. I also note that *you never actually attach the lead that you
created*; is this really what you want?

[1]: https://gitlab.kwant-project.org/kwant/kwant/issues/222

Re: [Kwant] About MoS2 ribbon

[Joseph Weston]

Hi,


> When I use kwant to calculate the conductance for MoS2 ribbon, zigzag
> edge is right, however, for armchair edge ribbon, errors occure as
> follows,
>
> LinAlgError: QZ iteration failed to converge in zgges.
>
> I am very comfused about this. Thanks!


Thanks for the report, other people have already reported a similar problem:

    https://gitlab.kwant-project.org/kwant/kwant/issues/176

Can you post a minimal example that demonstrates the error?

Happy Kwanting,

Joe

Re: [Kwant] Voltage

[Joseph Weston]

Hi Rohit,


> Is this code correct for the voltage applied in the 3rd lead?


I briefly looked at your code and you seem to be calculating the
scattering matrix for a stadium with 3 leads attached, except two of the
leads are 1 site wide and overlap. I am almost 100% sure this is not
what you intended, given the surrounding code.

If you provide a bit more context about what you're trying to
accomplish, any maybe ask a more focused question you're more likely to
get a useful response.

Happy Kwanting,

Joe

Re: [Kwant] regarding error message

[Joseph Weston]

Dear Shyam,

>
> File "C:\Users\Shyam Lochan Bora\Desktop\new 1.py", line 5, in 
>     from wraparound import wraparound
> ModuleNotFoundError: No module named 'wraparound'
>
> Help me in this regard

'wraparound' is now included with Kwant, as of version 1.3.
You can change the line:

    from wraparound import wraparound

to

    from kwant.wraparound import wraparound

to get wraparound working.

Happy Kwanting,

Joe

Re: [Kwant] Calculating valley polarized conductances using Kwant

[Joseph Weston]

Hi Kevin,

Thanks for the well-posed and focused question!

> I would like to calculate the valley polarized conductance for a graphene 
> nanoribbon with 2 leads (let’s say in the x-direction). My current approach 
> is this:
>
> 1) Get the S-matrix for a certain energy
> 2) Find the indices of the modes with positive velocities
> 3) Separate these modes by positive and negative momenta
> 4) Get the sub matrix from the S-matrix
> 5) Sum the transmission coefficients for positive and negative momenta

This is the right idea.

> 2) Find the indices of the modes with positive velocities

We only care about selecting (by their valley) the correct *incoming*
modes. In Kwant's convention these are the modes with *negative*, not
positive, velocity.

This is in an obscure part of the documentation [1], where we say

> positive velocity and momentum directions are defined with respect to
the translational symmetry direction of the [lead]

When leads are attached to a scattering region, their symmetry direction
points *away* from the scattering region "towards infinity", which means
that the incoming modes have negative velocity in Kwant's convention.


> The code below illustrates these 5 steps:
>
> import progressbar
> def valleyPolarizability(syst, energies, lead_start=0, lead_end=1):
> KPs = []
> Ks = []
> bar = progressbar.ProgressBar()
> for energy in bar(energies):
> smatrix = kwant.smatrix(syst, energy, params={'phi': phi})
> positives = np.where(smatrix.lead_info[lead_start].velocities >= 0)[0]
> momentas = smatrix.lead_info[lead_start].momenta[positives]
> K_prime_indices = np.where(momentas < 0)[0]
> K_indices = np.where(momentas >= 0)[0]
> submatrix = smatrix.submatrix(lead_end, lead_start)
> K_prime_T = np.absolute(np.sum(submatrix[:, K_prime_indices]))**2
> K_T = np.absolute(np.sum(submatrix[:, K_indices]))**2
> KPs.append(K_prime_T)
> Ks.append(K_T)
> return KPs, Ks

This is almost correct except for a couple of things. Firstly we should
select the modes that have negative velocity, as already discussed.
Secondly:

> K_prime_T = np.absolute(np.sum(submatrix[:, K_prime_indices]))**2

Here you're summing the scattering *amplitudes*, when you really want to
sum the scattering *probabilities*:

>  K_prime_T = np.sum(np.abs(submatrix[:, K_prime_indices]))**2

i.e. you need to take the absolute square before summing.


Happy Kwanting,

Joe


[1]:
https://kwant-project.org/doc/1.3/reference/generated/kwant.physics.PropagatingModes

Re: [Kwant] Onsite and hopping values in QHE

[Joseph Weston]

Hi


> sym = kwant.TranslationalSymmetry(
>     lattice.vec((1, 0, 0)),
>     lattice.vec((0, 1, 0)),
> lattice.vec((0, 0, 1))
> )
>
> []
>
> #Hall bar
> def onsite(site, B):
>   (x, y, z) = site.pos
>   return stored_model[site]
>
> def hopping_Ax(site1, site2, B):
>   x1, y1, z1 = site1.pos
>   x2, y2, z2 = site2.pos
>   return stored_model[site1,site2] * np.exp(-0.5j * B * (x1 + x2)
> * (y1 - y2))
>
> kwant_model[lattice.shape(shape, (0, 0, 15))] = onsite
> kwant_model[lattice.neighbors()] = hopping_Ax
>
> kwant_sys = wraparound.wraparound(kwant_model).finalized()
>
>  B = 0.02
> # Obtain the Hamiltonian as a dense matrix
> ham_mat = kwant_sys.hamiltonian_submatrix(args=[B], sparse=True)

The docstring for 'wraparound' [1] says that the wrapped around system will
have additional parameters, corresponding to the momenta. As you have a 3D
system and have wrapped around all 3 dimensions, the parameters will be
'k_x', 'k_y', and 'k_z'.

As you have only supplied a single argument (B) to your system, an error
is raised.
You will need to provide all the arguments to your system:

    kwant_sys.hamiltonian_submatrix(params=dict(k_x=0, k_y=0, k_z=0,
B=B), sparse=True)

Happy Kwanting,

Joe


[1]:
https://kwant-project.org/doc/1/reference/generated/kwant.wraparound.wraparound#kwant.wraparound.wraparound

Re: [Kwant] Band Structure in an Electric Field (Lead vs System)

[Joseph Weston]

Hi,


> I want to plot the band structure for my system, i.e. a quantum dot. I
> am still confused as to what code I need to write to plot that, or
> whether sys.leads[0] is actually correct for what I need.

Calculating a band structure for a quantum dot is not mathematically
well-defined.
A quantum dot does not have translational symmetry, so momentum is not a
good quantum number. You cannot label the eigenstates of your
Hamiltonian with a unique momentum 'k', so plotting the associated
eigenvalues as a function of 'k' (i.e. the band structure) is not possible.

Even if you attach leads to your quantum dot you'll still have the same
problem. While the leads *individually* have translational symmetry, the
combination of leads and quantum dot does not, because the quantum dot
breaks the translational symmetry.


Given this, it is difficult for us to interpret what you mean when you
say you want to plot "the band structure for your system". If you
provide more context for what your end goal is, we may be able to help more.

Happy Kwanting,

Joe

Re: [Kwant] Defining orbitals in a 3D structure with 3 basis atoms.

[Joseph Weston]

Hi again,

>
> But im still confused about some definitions
>
> 1) Alexandre, thank you for you suggestions, but honestly Im not
> pretty sure how to define a "superatom" that reproduces the structure
> im working on. How should i define the sites for each atom in the lattice?

He means that instead of defining a lattice with 3 sites per unit cell
with 3, 6, and 3 orbitals respectively, you should define a lattice with
a *single* site per unit cell, which has 3 + 6 + 3 = 12 orbitals. Of
course, plotting this thing is not so useful, and because you
essentially lose the information about where the "real" atoms are in
space (because you only have a single site representing 3 physical
atoms) it becomes more complicated to do things like define the action
of a position-dependent potential.

I would recommend keeping the physical atoms as separate sites (as you
are doing now).

>
> 2) Joseph, thank you for the detailed explanation, i rewrote that part
> of the code following your instructions and i think is better now, but
> i cant figure a way to describe the interactions that occur "inside"
> the same atom, for example, Re has 6 orbitals, how can i express the
> hoppings that have place in the same site?

If a site on one of your lattices represents a Re atom with 6 orbitals,
then what Kwant calls the "onsite" Hamiltonian terms must be 6x6
matrices: the diagonal contains the actual onsite values, and the
off-diagonal terms are hoppings between orbitals of the *same* Re atom
(this is why Kwant still calls this an "onsite").

Similarly, hoppings between *different* Re atoms must be 6x6 matrices,
and hoppings between Re and C (or N) atoms must be 6x3 matrices, to
match the number of orbitals on each end of the hopping.

In your code example there is a section where you set the onsite and
hoppings for the 3 atoms inside the unit cell. This is nearly correct,
except that  you have  12x12 matrices for each of the onsite/hopping
values, when you really need matrices with shapes that match the number
of orbitals on the sites.

On the other hand, you then use 'lat.shape' to assign an onsite value of
'4' to all sites, so I think that you are confused about what exactly
you are doing. You will have to use 'C.shape', 'Re.shape', and 'N.shape'
(C, Re and N being your sublattices) to assign the onsite matrices of
the appropriate shape to sites from each of the sublattices in turn,
because sites from each sublattice have different numbers of orbitals.

> P.D. The code so far is giving me an error in the leads apparently
>
> /home/sergio/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:197: 
> RuntimeWarning: When finalizing lead 0: Infinite system with disconnected 
> cells.
> /home/sergio/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:197: 
> RuntimeWarning: When finalizing lead 1: Infinite system with disconnected 
> cells.

As the warning illustrates, the unit cells of your leads are not
connected by hoppings, so all elements of the scattering matrix will be
0. Looking at your code I see that you indeed don't add any hoppings to
your leads; you only populate a bunch of sites using 'lat.shape'.

Happy Kwanting,

Joe

Re: [Kwant] Onsite and hopping values in QHE

[Joseph Weston]

HI,


>
> I am trying to reproduce the QHE bar device calculations
>
> http://nbviewer.jupyter.org/github/topocm/topocm_content/blob/master/w3_pump_QHE/Laughlinargument.ipynb
>
>
> However, I am having difficulty with the onsite and hopping values as
> I am using TBModels
>
> https://github.com/Z2PackDev/TBmodels
>
> and therefore the values are not the same for all sites.

I am unsure of the output format of TBModels. Glancing at their
documentation I can't see any way of getting the raw matrix elements
from a "Model". Looking into the source code it seems that they are
stored in a format that will be annoying to query manually.

> The problem is that after the scalar values have been passed on to the
> Builder, these will then be overriden when using functions for the
> position-dependent Hamiltonian.
>
> I tried making a deep copy of the model, which would then serve as a
> data structure where to draw the scalar values from when calling the
> Onsite and Hopping functions. The problem with this is
>
> a) I am getting a KeyError which I am still not sure why it happens
> b) It does not feel very memory-efficient.
>
> I am wondering if there is another way that this could be done on the
> Kwant side that I have not been able to think of.

I can't think of a better way to accomplish this. Using TBModels'
internal datastructures directly is IMO too much hassle unless your
systems are really big.

Did you check that the symmetry of your Kwant system and the TBModels
Model are compatible?

Joe

Re: [Kwant] Josephson current in a lattice between two superconductors

[Joseph Weston]

HI Ville,


> I'm trying to use Kwant to investigate behavior of a lattice between
> two superconductors, namely trying to figure out the Josephson
> current. The task is basically to diagonalize the related
> Bogoliubov-de Gennes Hamiltonian and use the eigenvectors and
> eigenvalues to calculate the current in the lattice. I basically
> follow the example 2.6. in Kwant documentation to build up the
> Hamiltonian but let the both leads be superconductors and establish
> the system in between. I introduce the phase difference between the
> leads by multiplying the Delta parameter of the other lead by e^(i
> phi). However, I'm a bit stuck there as I don't know how to utilize
> the obtained eigenvectors to make wave functions which I could input
> to t current operator in Kwant.
>

To get the correct Josephson current you'll need to include the
contribution from the Andreev bound states that form in the Junction. As
these are true bound states you won't be able to get them by using
'kwant.wave_function', as the latter only calculates the *scattering*
states.

Some of the Kwant Authors recently published a paper that describes how
to calculate bound states for general scattering systems [1] (i.e.
systems with leads). We're working to implement this in Kwant, but it's
not there yet.

As far as I am aware the following are possible approaches to calculate
the bound state spectrum are:

+ truncate the leads after some (large) length and diagonalize the
resulting system. You'll have to identify "by eye" which states
correspond to the true bound states of the infinite system you are
approximating.  The leads will need to be truncated far enough away from
the scattering region that the bound state wavefunction is nearly vanishing.

+ Replace the leads with a self-energy term and diagonalize the
resulting matrix. I can't remember what guarantees one gets about using
this approach, as you make your Hamiltonian non-Hermitian by doing this.

+ Construct your system as 3 parts: an SN interface, an N-N interface
(include any disorder or whatever here), and an NS interface. Calculate
the scattering matrices of the 3 parts separately, and then the bound
states can be found by combining the scattering matrices in a particular
way and solving for zero determinant (I can't remember the details, but
it's in section 2 of [2], and references therein).


Basically everyone in our group in Delft calculates Josephson currents
as a matter of course, so they're probably in a better position to tell
you which way to go. I thought I'd give my 2 cents anyway.

Happy Kwanting!

Joe

[1]: https://arxiv.org/pdf/1711.08250.pdf
[2]: https://arxiv.org/pdf/cond-mat/0406127.pdf

Re: [Kwant] Coupling Matrix between Lead and Conductor

[Joseph Weston]

Hi,


>
> I have a matrix giving me the complex wave function at each site. How
> do I know which mode (m/n) do the amplitudes correspond to?

Where did you get this matrix from? Is this the wavefunction of a mode
in a lead, from 'kwant.modes', or of a scattering wavefunction, from
'kwant.wave_function'? Bear in mind that a mode wavefunction and a
scattering wavefunction are defined in different places: a mode
wavefunction is defined in the lead only (i.e. over a single unit cell,
and you use Bloch's theorem to get the wavefunction in neighboring unit
cells), whereas the scattering wavefunction may only be queried in the
scattering region itself (even though it extends infinitely far into the
leads, Kwant does not allow you to query the amplitude of a scattering
wavefunction on sites the leads).

> Also, I am not aware of the method to get wavenumbers k_m. Could you
> kindly help me out with that as well?
>
You can get the wavenumbers by calling 'kwant.modes' on a lead and
querying the propagating modes. This contains the wavefunctions for the
modes of the provided lead, as well as the wavenumbers and velocities.

Happy Kwanting,

Joe

Re: [Kwant] Defining orbitals in a 3D structure with 3 basis atoms.

[Joseph Weston]

Hi all,


Sergio, if you want to create a Polyatomic lattice you have to pass all
the basis vectors and orbital numbers in a single call to
'kwant.lattice.general', otherwise there's no way Kwant can know that
you want the 3 sublattices to be considered as a single Polyatomic lattice!

    prim_vecs = [...]  # Bravais lattice vectors

    basis [ [...], [...], [...]]  # position vectors of the 3 basis atoms

    norbs = [3, 6, 3]

    lat = general(prim_vecs, basis=basis, norbs=norbs)

    a, b, c = lat.sublattices

    ...

    syst[a(0, 0, 0), b(0, 0, 0)] =  [...] # a-b hopping in same unit cell
    syst[a(0, 0, 0), c(1, 0, 0)] =  [...] # a-c hopping in between unit
cells

This is covered briefly in the tutorial [1]; do you think we can
improved that section of the tutorial in some way?

Happy Kwanting,

Joe


[1]: https://kwant-project.org/doc/1/tutorial/graphene



On 06/07/2018 12:05 PM, alexandre.berna...@u-psud.fr wrote:
> Hello,
>
> First of all, I am new to kwant, so I hope I will give you relevant
> insights.
>
> It seems to me that:
>
> 1) There are some mistakes in the syntax you are using, so you should
> check again the examples available in the doc :
> https://kwant-project.org/doc/1/
>
> 2) Since you consider different numbers of orbitals for your atoms, I
> think you should use only one "superatom" with 12 orbitals per unit cell.
> Then Kwant allows to define 12x12 matrices both for the on-site terms
> (your matrix elements within the superatom) and for the hopping terms
> (your matrix elements between superatoms).
>
> I hope this is a bit helpfull.
>
> Best regards,
>
> Alexandre BERNARD
>
> 
> *De: *"Sergio Castillo Robles" 
> *À: *kwant-discuss@kwant-project.org
> *Envoyé: *Mercredi 6 Juin 2018 03:54:21
> *Objet: *[Kwant] Defining orbitals in a 3D structure with 3 basis atoms.
>
> Hello everyone, i would appreciate any suggestion you could give me to
> solve this.
>
> Well, im trying to create a structure with 3 basis atoms in the unit
> cell, each atom has a different number of orbitals a=3, b=6 and c=3. I
> have tried with polyatomic module and creating one lattice for each
> atom with no success
>
> The thing is that i cant figure it out how to create a lattice with
> different orbitals in each site. I have the on-site and hopping
> energies matrix so i need to be able to introduce this values.
>
> Im pretty new using kwant so any advice is appreciated.
>
> Here is the code i've working on (ignore the matrix elements for the
> moment):
>
> import kwant
> import tinyarray
> import numpy
> from matplotlib import pyplot
>
> VecPrim = [(3.176, 0, 0), (1.588, 2.7504, 0), (0, 0, 17.49)]
> base =  [(0, 0, 0), (1.588, 0.9168, 0.8745), (0, 0, 1.749)]
>
>    ### This is the part that i cant get it right, i have tried
> kwant.lattice.polyatomic
>    ### and defining 3 lattices one for each atom, but no success
>    
> a = kwant.lattice.general(prim_vecs = VecPrim, basis = (0, 0, 0))
> Cpx, Cpy, Cpz, = a.sublattices
> b = kwant.lattice.general(prim_vecs = VecPrim, basis = (1.588, 0.9168,
> 0.8745))
> xy, xz, yz, x2y2, z2, s = b.sublattices
> c = kwant.lattice.general(prim_vecs = VecPrim, basis = (0, 0, 1.749))
> Npx, Npy, Npz = c.sublattices
>
>     ### Ignore from here ###
>
> def make_cuboid(t=1.0, a=15, b=10, c=5):
>     def cuboid_shape(pos):
>     x, y, z = pos
>     return 0 <= x < a and 0 <= y < b and 0 <= z < c 
>     syst = kwant.Builder()
>     syst[lat.shape(cuboid_shape, (0, 0, 0))] = 4 * t
>    
>     syst[lat(0, 0, 0)] = numpy.array([[-1.42, -0.014064-0.000409j,
> 0.026908+0.003422j, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [-0.014064+0.000409j, -1.784936,
> -0.031384+0.021565j, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0.026908-0.003422j,
> -0.031384-0.021565j, 0.710443, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> 0],  
>  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> 0]])
>    
>     syst[lat(0, 1, 0)] = numpy.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0,
> 0, 0, 0],
>  [0, 0, 0, 0, 0, 0, 0, 0, 0,
> 0, 0, 0],
>  [0, 0, 0, 7.168134,
> -0.005189-0.047376j, 

Re: [Kwant] Changing the potential in a continuum Hamiltonian

[Joseph Weston]

Hi,

You can do what you want by declaring extra parameters to the 'V'
function when you make your Hamiltonian:


    H = 'k_x**2 + k_y**2 + V(x, y, Vgate)'

    template = kwant.continuum.discretize(H)

    syst = kwant.Builder()

    syst.fill(template, ...)  

    ...

    def gate(x, y, Vgate):

    return ...

    ...

    syst = syst.finalized()

    kwant.smatrix(syst, params=dict(V=gate, Vgate=3))


Note how 'Vgate' appears in the original Hamiltonian string, as a
parameter to 'gate', and as an argument when passing parameters to the
system when calling 'smatrix'.

Happy Kwanting,

Joe


On 06/04/2018 05:15 PM, jonas.wiedenm...@physik.uni-wuerzburg.de wrote:
> Hi everybody,
>
> I hope my problem is rather simple to solve/explain.
> I would like to create a rectangular system using the
> function kwant.continuum.discretize, apply a local potential, and plot
> the conductance of the system as a function of this potential.
> My Problem is: How to pass the function potential to the system. 
>
> So far I used a similar version of the code from tutorial 2.10 
>
> /
> def system(L=70, W=40):
>     
>     hamiltonian = "k_x**2 + k_y**2 + V(x, y)"
>     template = kwant.continuum.discretize(hamiltonian) 
>    
>     def shape(site):
>         (x, y) = site.pos
>         return (0 <= y < W and 0 <= x < L)
>     /
>     
> /
>     def lead_shape(site):
>         (x, y) = site.pos
>         return (0 <= y < W)
>     syst = kwant.Builder()
>     syst.fill(template, shape, (0, 0));
>     lead = kwant.Builder(kwant.TranslationalSymmetry([-a, 0]))
>     lead.fill(template, lead_shape, (0, 0))
>     syst.attach_lead(lead)
>     syst.attach_lead(lead.reversed())
>     /
> kwant.plotter.map(syst, lambda s: V(s, 0.5));
> /
>     /
>    
> /
>     return syst
> /
>
> Then, I defined a potential which returns for a given position x,y a
> number (This is taken from a kwant tutorial named: Transport through a
> barrier)
>
> /
> def rectangular_gate_pot(distance, left, right, bottom, top):
>     
>     d, l, r, b, t = distance, left, right, bottom, top
>
>     def g(u, v):
>         return atan2(u * v, d * sqrt(u**2 + v**2 + d**2)) / (2 * pi)
>
>     def func(x, y, voltage):
>         return voltage * (g(x-l, y-b) + g(x-l, t-y) +
>                           g(r-x, y-b) + g(r-x, t-y))
>
>     return func
>
> gate = rectangular_gate_pot(10, 30, 40, -10, 50)
>
> /
>
> /
> def potential(site, V):
> /
>     x, y = site.pos
>     return gate(x, y, V) 
>
> Next, I would like to calculate the conductance. I used again a
> snipped from the tutorial but I do not understand how I can properly
> address the potential
>
> /def plot_transmission(syst, energy, voltages):/
> /    trans = []/
> /    for voltage in voltages:/
> /        smatrix = kwant.smatrix(syst, energy, /dict(V=potential(*what
> do I have to put in here or does it not work this way??*)/)/
> /        trans.append(smatrix.transmission(1, 0))/
> /    pyplot.plot(params, trans)/
>
>

Re: [Kwant] seeking help

[Joseph Weston]

Hi,


> Dear Sir,
>   
>      I want to model a bi2se3 nanowire in kwant to investigate the
> transport properties.
>    
> Please help me in this regard.

Your question may have several interpretations:

1. How do I make a tight-binding model for Bi2Se3 that includes all
relevant physics for my application?
2. How do I write down a (known) tight-binding model in Kwant, in general?
3. How do I write down a (known) tight-binding model in Kwant,
specifically for my nanowire setup?

If your question is really (1) then I am afraid that this mailing list
is probably not the right place
to be asking, as this mailing list is about Kwant-specific questions.

If your question is really (2) then I recommend going through the Kwant
tutorial [1].

If your question is really (3) then I recommend going through the Kwant
tutorial [1]. If there are
some intricacies in your model that you don't know how to express in
Kwant, then by all means
post back to this mailing list with your specific question.

Happy Kwanting,

Joe


[1]: https://kwant-project.org/doc/1/tutorial/

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

[Joseph Weston]

Hi again,

> Hello Joseph,
>
> One problem that I see is that I can only get neighbors() from
> Builders and lattices (using site.family).
>
> Using something like this:
>
> def first_neighbors(to_site, from_site):
>   return to_site in from_site.family.neighbors(1)
>
> is a good idea, but am I really comparing two sites?
>
> 
> for i in x.family.neighbors(1):
>     print(i)
>
> HoppingKind((1, 0, 0), )
> HoppingKind((0, 1, 0), )
> ---
>

Ah, sorry I was confused.


> Also, since this is very expensive computationally, is there any other
> way to find the hopping distance between two sites?

Sure! If the sites are from the same lattice then you can just take the
difference of the "tags" of the sites:

    delta = to_site.tag - from_site.tag

'delta' will be an integer array: the position of the site in the basis
of lattice vectors. If the sites are from different
lattices then you probably can't do better than comparing the realspace
positions:

    realspace_delta = to_site.pos - from_site.pos

Happy Kwanting,

Joe

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

[Joseph Weston]

Hi,

>
> Is it also possible to help me with this function to match first
> nearest neighbor sites for the where argument in kwant.operator.Current:
>
> def first_neighbors(site_to, site_from):
> for i in site_from.family.neighbors(1):
>   if(??):
>     return True
>   else
>     return False
>
> [...]
>
> J = kwant.operator.Current(kwant_sys, where=first_neighbors, sum=True]

Probably the following:

    def first_neighbors(to_site, from_site):
        return to_site in from_site.neighbors(1)

'neighbors' returns (per the documentation [1]) a list of sites and we
can just use
the 'in' operator to check for the requested site.


Happy Kwanting,

Joe


[1]:
https://kwant-project.org/doc/1.0/reference/generated/kwant.lattice.Polyatomic#kwant.lattice.Polyatomic.neighbors

Re: [Kwant] Four leads or six leads in hall bar-How to use " add_site_family" exactly in Kwant

[Joseph Weston]

Hi Bill,


> Dear Joe   
>     "add_site_family" did work in graphene zigzag nanoribbon.When I
> use "add_site_family",the central scattering region would not enlarge
> after connected with leads in zigzag nanoribbon.But when I use it in
> graphene armchair nanoribbon,the central scattering region will
> enlarge. If you have spare time,you can run my code,and you can see
> the picture that the central scattering is enlarge a little in up and
> downd leads,left and right lead is right.
>       Thank you
>        Bill Yang
>


Ah, I understand now. In your original email you talked about a
"triangle of sites", which I did not see in the images you posted, so I
was unsure what you were talking about.

I don't think you be able to get rid of the "extra added bit" with a
judicious choice of 'other_vectors'; I think you'll have to define the
lattice differently in order to do what you want.

On the other hand, is this extra added region such a problem? There will
be no physical consequences, as the extra added piece has the same
Hamiltonian as the leads...

Happy Kwanting,

Joe

Re: [Kwant] Four leads or six leads in hall bar-How to use " add_site_family" exactly in Kwant

[Joseph Weston]

Dear Bill,


> This is David's question.Link is 
> https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2014-October/000169.html
> .

In a response to that question Michael Wimmer posted example code [1] that uses 
'add_site_family'. Did you try doing as he suggests in that post?

>   I use "add_site_family",but it didn't work.

What, exactly, did not work? Did it raise an exception? Did it not raise
an exception, but the result was not what you expected?


Happy Kwanting,


Joe


[1]:
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2014-October/000174.html

[Kwant] Kwant extension directory

[Joseph Weston]

Hello everybody!


We recently added a new section to the Kwant website [1] that contains
links to pieces of code that "extend" Kwant in some cool way. If you've
used Kwant in a way that you feel may be useful to the wider community
make sure to get in touch so that we can add it to the website. The best
way to do that is to make a post to the Kwant development mailing list
[2] with a link to your code and a few sentences describing what it
enables you to do.


Happy Kwanting,


Joe


[1]: https://kwant-project.org/extensions

[2]: mailto:kwant-de...@kwant-project.org

Re: [Kwant] NNN hopping with periodic BC

[Joseph Weston]

Sorry I only just saw this email; please remember to "reply All" so that
the email gets sent to the mailing list too!

> Thanks. But except the PBC imposed on transverse direction, I need
> attach the lead on longitude direction. So I want the transverse size
> of sample to be finite.
>

Ah I see. In this case your could consider using 'kwant.wraparound' [1].
In your case you would want to set the 'keep' paramter to '0' to keep
the 0th translational symmetry in the wrapped around system. Using
wraparound in this way will produce a Builder with a single
(longitudinal) translational symmetry that you may use as a lead. The
produced Builder will have an (additional) parameter 'k_y'; the
transverse momentum. Don't be confused by the name; it will also work
when the "wrapped" symmetry vector was not in the "y" direction!

There were already several people on the mailing list who used
wraparound, so you could also search the archives. If it's still not
clear after this post back to this thread.

Happy Kwanting,

Joe

[1]:
https://kwant-project.org/doc/1/reference/generated/kwant.wraparound.wraparound#kwant.wraparound.wraparound

Re: [Kwant] NNN hopping with periodic BC

[Joseph Weston]

Hi,

> I want to create a graphene-like system with periodic BC at transverse
> direction.  I’m following
>
> the discussion at
> https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00036.html.
>
>  
>

If you only want 1 direction to have periodic boundary conditions (and
finite in the other direction)
why not just create a system with a 1D translational symmetry? This will
make your life a lot easier.

Happy Kwanting,

Joe

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

[Joseph Weston]

Hi,


>
> So, for example, I can print the on-site energies or the hopping
> integrals of the sites I want using:
>
> ..
>
> However, I am still unsure about the 'where' argument
>
> -
>
>
 J = kwant.operator.Current(kwant_sys,
 where=kwant_sys[lattice.sublattices[4](0,0,-1),
 lattice.sublattices[6](0,0,-1)])
> Traceback (most recent call last):
>   File "", line 1, in 
>   File "kwant/operator.pyx", line 883, in kwant.operator.Current.__init__
>   File "kwant/operator.pyx", line 217, in
> kwant.operator._normalize_hopping_where
>   File "kwant/operator.pyx", line 218, in genexpr
> ValueError: need more than 1 value to unpack
 J = kwant.operator.Current(kwant_sys,
 where=(lattice.sublattices[4](0,0,-1),
 lattice.sublattices[6](0,0,-1)))
> Traceback (most recent call last):
>   File "", line 1, in 
>   File "kwant/operator.pyx", line 883, in kwant.operator.Current.__init__
>   File "kwant/operator.pyx", line 223, in genexpr
>   File "kwant/operator.pyx", line 223, in genexpr
> AttributeError: 'Builder' object has no attribute 'graph'
>

There are 2 problems here:

1. 'where' expects a *sequence* of places where you would like to
evaluate the current. In your code snippet you pass the value function
associated with some hopping. Using your example this would be something
like this:

    lat1 = lattice.sublattices[4]
    lat2 = lattice.sublattices[6]
    J = kwant.operator.Current(syst, where=[(lat1(0, 0, -1), lat2(0, 0,
-1))])
  

2. 'Current' expects a *finalized* system, but you pass it an
*unfinalized* system (Kwant calls this a "Builder").


Does that make sense? Can you also please point me to where you looked
in the documentation and explain which bits were confusing? This will
help us improve the documentation!

Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] modes and current density

[Joseph Weston]

Hi,



> 
> In my case, I have a system with more than one mode. I have found some
> lines that consider all the modes (at the end of this page :
> https://kwant-project.org/doc/dev/tutorial/operators)
> 2.7. Computing local quantities: densities and currents ...
> 
> kwant-project.org
> 2.7. Computing local quantities: densities and currents¶ In the
> previous tutorials we have mainly concentrated on calculating global
> properties such as conductance ...
>
>  
> 
> wf_left=wf(0)
> J_m_left=sum(J_m_bound(p)forpinwf_left)
> J_z_left=sum(J_z_bound(p)forpinwf_left)
> 
> But this is equivalent to calculate the expression
>
>
> which is simply the sum over the current coming from each mode.
> However, I would have rather use the following expression to find a
> more general expression of the current density
>  
>
>
> What is therefore the good way to find J? I have not found more
> informations in the Kwant documentation.
>
> Best regards,
>
> Nicolas

If you want to calculate the total current at a given energy you really
should sum the contributions from each mode independently!

You say that the "more general" expression for the current density is:


but I'm not sure what you mean by that. In the above expression you are
making some weird kind of mix of different modes. This is certainly not
the expectation value of the current operator; are you maybe trying to
calculate some kind of noise properties?

Regards,

Joe

Re: [Kwant] Building tight binding Hamiltonian with multi atom with multi-orbital basis

[Joseph Weston]

Hi,


> Thank you for your suggestion. Understanding the way you said, I made
> the tutorial code for band_structure.py with 4 x 4 Hamiltonian
> (attached to this mail). I have defined tau_x  and tau_z as 4 x 4
> Hamiltonian as done in superconductivity.py tutorial. I do not know
> whether that is the right way to make multiorbital code here or not.
> So, please tell me looking at the attached code.

Your script correctly constructs a multi-orbital system. I am not sure
what this system corresponds to physically, but this is a question that
only you can answer because it is your model.

> But the problem is how to specify which two orbitals belong to one
> site and other two orbitals belong to another site (suppose there are
> two sites)?

I am not sure what you mean. When you write a line like

    syst[lat(i, j)] = tau_z

you are saying that the Hamiltonian on site 'lat(i, j)' is the (matrix)
value 'tau_z'. This implies that you have several orbitals (4 in your
case) on site 'lat(i, j)'. There is no longer the concept of your 4
orbitals being "built" from others (e.g. s & p orbitals "outer
producted" with spin up & spin down); you have to keep track of what
your orbitals mean.


> Also if there is any way to give orbital symmetry specifically, like
> px, py or dxy etc.

You will have to specify this symmetry yourself by making sure that your
matrix elements satisfy the appropriate symmetry.


> Even If eigenstate of the Hamiltonian can be plotted in real-space ?
You can plot, for example, the square magnitude of a wavefunction using
'kwant.plotter.map'; is this what you mean?

I'm not sure if I clarified anything for you, but if you still have
doubts don't hesitate to post back in this thread.

Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Compute conductance for different system lengths

[Joseph Weston]

Hi Jan,


>
> I'm looking for a workaround for a problem I currently face: when
> computing the two-terminal conductance for a system of length L, is
> there a way to get the conductance for the same system of length
> 1,2,..L-1 on the fly (while keeping everything else the same, i.e.,
> same width, same disorder configuration etc.)?

Kwant does not support this. You would need to construct a new system
for each of these cases, which would incur the relatively high cost of
system construction and finalization for every value of L. A possible
workaround would be to construct the scattering system of length L, and
then to add a parameter to your onsite/hopping functions that you can
tune, such that your "effective" scattering region is whatever length
you like. This would be faster, as you would only have to construct and
finalize your system once.

> As far as I understand how the scattering matrix calculation works
> internally, it shouldn't take much longer to compute these
> intermediate values than just getting the final conductance.

I'm not sure what you mean. The default scattering solver does not use
the recursive Green's function method or anything. We (more or less)
solve a linear system in the basis of modes in the lead, and local
degrees of freedom in the scattering region, so that the total solution
vector contains scattering matrix components in the "lead" part, and the
scattering wavefunction in the "scattering region" part.  We set up this
linear system and then pass it off to a sparse linear solver (MUMPS by
default). Given this, it is not immediately obvious to me how we would
compute these "intermediate" values.

Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Importing Hamiltonian from file

[Joseph Weston]

Dear Marta,

> 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?

I am not sure what you are exactly asking. You would need to read in the
Hamiltonian from the file and assign the onsites and values to a
Builder, as you would with any other Kwant system. Am I missing something?

Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Mistake on the docstring of `PropagatingModes` ?

[Joseph Weston]

Hi,


The documentation was corrected, but the new documentation version is
not in the "stable" documentation yet because we did not make a new
release. The updated documentation is available in the "dev" docs:

https://kwant-project.org/doc/dev/reference/generated/kwant.physics.PropagatingModes


Happy Kwanting,


Joe


On 10/16/2017 11:23 AM, George Datseris wrote:
> Dear All,
>
> (thanks so much for kwant ;) )
>
> I am wondering whether there is a mistake/typo in the docstring of the
> PropagatingModes object.
>
> The reason is the following: In a simple lead (specifically graphene
> zigzag nanoribbon) at energy that excites 6 modes (only 3 with
> positive derivative dE/dk),
> when getting the `PropagatingModes` (either through physics.modes or
> by the Smatrix), and I get the momenta and velocities, it looks like
> this:
>
> In [35]: modes
> Out[35]: 
>
> In [36]: modes.momenta
> Out[36]:
> array([-2.12874239,  2.0566145 ,  2.06992722,  2.12874239, -2.0566145 ,
>    -2.06992722])
>
> In [37]: modes.velocities
> Out[37]:
> array([-1.91098658, -2.24527103, -1.66541483,  1.91098658, 2.24527103,
>     1.66541483])
>
>
> Now, the docstring of Propagating modes states:
>
>> The first half of the
>> modes have negative velocity, the second half have positive velocity.
>> The
>> modes with negative velocity are ordered from larger to lower
>> momenta, the
>> modes with positive velocity vice versa.
> However, this is not what I see. The first half, with negative
> velocity, is clearly sorted from smaller to bigger momenta, while the
> second half (with positive velocity) is sorted from bigger to slower
> momenta. This is the exact opposite of the docstring.
>
> Thus I am writing this mail to see if I have misunderstood something
> or if this is simply a typo mistake.
>
> Best,
> George
>




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Band structure with a periodicity larger than the lattice constant

[Joseph Weston]

Dear Michael,


> sorry for bothering again on this problem, but implementing what
> written in the previous email I am getting some unexpected results.
> Using /d=4*a/ and not including anything with a periodicity different
> from the one of the honeycomb lattice the band structure quadruple and
> each duplicate shifts of \pi/2 as in the figure below. What am I missing?

It seems to me that you just get band folding, as expected. When you
increase the symmetry period you reduce the size of the Brillouin zone,
but you increase the number of bands (because there are now more sites
in the unit cell).

What did you expect?

Happy Kwanting,

Joe


signature.asc
Description: OpenPGP digital signature

Re: [Kwant] definition of tag

[Joseph Weston]

Dear KS,


> I want to know more about Hoppingkind and tag. Is the tag of a site
> the order pair (i,j ) in  lat(i,j) which is the site considered? This
> is for two dimension, for three dimension it should be (i,j,k) of lat
> (I,j,k).
>

In the case you describe, the tag is indeed the ordered pair (i, j) and
(i, j, k).

In Kwant sites are uniquely identified by their "family" and a "tag". In
practice the "family" corresponds to the lattice that the site belongs
to, and the tag is an integer vector; the position of the site in the
basis of lattice vectors. Lattice objects can be called like functions,
with the tag as argument, and return the appropriate Site object.

Kwant actually defines "tag" in a more general way (because site
families in general need not be Bravais lattices ), but in practice most
users only use sites belonging to lattices (Bravais lattices and
collections of such; polyatomic lattices).

Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] Band structure with a periodicity larger than the lattice constant

[Joseph Weston]

Dear Michael,


> I would like to compute the band structure of a honeycomb lattice with
> a periodic field of periodicity d, with /d>a /(/a/ the lattice
> constant). What I have thought is to impose a translation symmetry to
> the leads with a period /d/ instead of /a/, as written below. Can
> kwant do this

Kwant does, of course, allow you to do this. You only need to make sure
that the translational symmetry you define is commensurate with the
lattices of any sites you add to the system (i.e. that the symmetry
operations are also valid lattice translations).


> is this the correct way to do it?
>
>     sym0 = kwant.TranslationalSymmetry((-d, 0.))
>     sym1 = kwant.TranslationalSymmetry((d, 0.))

This looks correct, although I am not sure why you are defining two
symmetries in opposite directions. If you want to add sites from the
honeycomb lattice to a system that has this symmetry you will need to
ensure that the symmetry is commensurate with the lattice (i.e. that
realspace vector (-d, 0) is an integer linear combination of the
honeycomb lattice vectors).


>
>     sym0.add_site_family(A, other_vectors=[(-d, 2*d)])
>     sym0.add_site_family(B, other_vectors=[(-d, 2*d)])
>
>     sym1.add_site_family(A, other_vectors=[(-d, 2*d)])
>     sym1.add_site_family(B, other_vectors=[(-d, 2*d)])

I am not entirely sure what you want to do here. "add_site_family" is a
low-level functionality of symmetries that is related to how Kwant will
choose the unit cell of the symmetry. It is usually used when defining a
lead that will later be attached to a finite scattering region, to avoid
"ears" getting added to the scattering region (see the images in
https://kwant-project.org/doc/1/tutorial/graphene for an example of what
I mean).

As all you want to compute is band structure (no attaching), I am unsure
as to why you would need this.

Happy Kwanting,

Joe


signature.asc
Description: OpenPGP digital signature

Re: [Kwant] A weird result of 1D quantum wire

[Joseph Weston]

Dear Liu,


>   I have met a strange problem in computing conductance of 1D
> quantum wire.
>   The setup is as follows:
>   The hopping of the quantum wire is given by the first kind of
> Bessel function, for example, J0(A).
>   The coupling of the lead is set as 1.
>   As we know in this case the transmission should be proportional
> to J0(A)^2. However, I get very different results for wire of even and
> odd sites numbers as shown in the attachment...
>   Is this a bug or do I miss something else?
>

This is probably not a bug in Kwant. I see in your script that you do
not construct your system in a very idiomic manner (you manually add
nearest-neighbor hoppings etc). Are you sure you did not make a mistake
during system construction? Did you plot the system in the two cases
(even/odd) and see that everything is as you expect?

Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature

Re: [Kwant] beyond effective mass (limiting kwant to a range of k)

[Joseph Weston]

Hi Leon,


Rafal's got it right, I'm afraid. The easiest way to solve this problem
is probably to modify your model so as to "push away" the modes with
large momentum so that they have energies outside of the relevant window
that you are concerned with.


Unfortunately you can't just "throw away" the modes at large momenta,
because this would make the scattering problem under-determined.
Although one could maybe find the "best" solution to this
under-determined problem using least squares or something, this is quite
awkward and I'm not even sure it's valid.

In order to properly "throw away" the unwanted modes, you would instead
have to make them evanescent. Doing this is essentially the same as
modifying your model so as to "push the modes away".


Happy Kwanting,


Joe




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] lead self energy

[Joseph Weston]

Hello Chung,


> 1-The lead's self energy is in principle a complex number. In kwant,
> If I want to add the self energy of a lead by hand, I can add the real
> part of it to the on-site energy of the site that lead is connected
> to. Where should I add the imaginary part of the self energy?
Kwant allows to to specify a lead by its self energy if you want. You
need only construct a "SelfEnergyLead" and manually attach it to your
Builder:

lat = kwant.lattice.chain(norbs=1)

syst = kwant.Builder()
syst[(lat(i) for i in range(4))] = 4
syst[lat.neighbors()] = -1

   # not sure what this corresponds to, but whatever
def self_energy(energy, args=(), *, params=None):
return -1j

# 'interface' is where we wish to "attach" the lead to the system
left_lead = kwant.builder.SelfEnergyLead(self_energy,
interface=[lat(0)])
right_lead = kwant.builder.SelfEnergyLead(self_energy,
interface=[lat(3)])

syst.leads.append(left_lead)  # tell the system about the lead
syst.leads.append(right_lead)

syst = syst.finalized()

In the above I made the assumption that you have a 1D system. In the
more general case 'self_energy' would need to return a matrix. Also the
self energy chosen above does not really correspond to anything in
particular (the self energy is energy dependent in general), it is just
meant to be illustrative.

You could also just put complex onsites on the interface sites in your
system directly, but then you would not be able to calculate transport
properties using Kwant (because you would not have told Kwant that there
were any leads!)


> 2- If I know the value of the level width function (\Gamma = 2πρt^2 =
> 1), and the lead is a 1D chain, can I calculate the value of the
> coupling parameter, t? Is there a way to obtain the DOS of the lead
> and then get t?
Sure, the DOS is just the reciprocal of dE/dk, the derivative of the
dispersion relation in the lead (up to a bunch or prefactors etc.)

Happy Kwanting,

Joe




signature.asc
Description: OpenPGP digital signature

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

[Joseph Weston]

Hi Mitchell,

> Sorry again for the reply, I've had a look at the documentation and it
> seems that relwidth and abswidth control the size of the arrays. What
> I'm stuck on though is that my system has   L = 80, W = 450 and H =52,
> I was expecting field to return something of the form (80,450,52,3), I
> must apologise again, I'm pretty new to all things python so my
> knowledge in the more complex areas is lacking.

The input to 'interpolate_current' is a current defined on the hoppings
of your system. The output is a current defined over a realspace grid.
First the current defined on the hoppings is convoluted with a bump
function to smooth it out, then the resulting continuous vector field is
sampled on a regular grid. The docstring explains how this is done. The
interpolation grid does not in general coincide with the "grid" of sites
that form your system. Typically there will be several interpolation
points between neighboring sites.

The 'box' returned by 'interpolate_current' allows you to calculate the
realspace position of the interpolation points.

Also the modifications to 'interpolate_current' so that it functions
correctly for 3D systems just landed on the most recent development
version of kwant [1].


Happy Kwanting,

Joe


[1]:
https://gitlab.kwant-project.org/kwant/kwant/commit/be364e7c64c316caf86869df397446399bcef1af




signature.asc
Description: OpenPGP digital signature

Re: [Kwant] continuous vs discrete leads

[Joseph Weston]

Dear Patrik,


> When I replace the 'if continuous' section with the following section
> it will complain.
>
> if continuous:
> def lead_shape(site):
> (x, y, z) = site.pos
> return (x == -0.7 and y == 0.0 and z <= -a)
> t00=0.0
> Leadham
> ="t00*sigma_0*k_x**2+t00*sigma_0*k_y**2-t0*sigma_0*k_z**2+(2*t0+e1)*sigma_0"
> template = kwant.continuum.discretize(Leadham, grid_spacing=a)
> dn_lead.fill(template, lead_shape, lat1(0, 0, -1))
> syst.attach_lead(dn_lead)

The error message I got when running the code was:

> ValueError: Lead to be attached contains no sites.

This error seems perfectly clear to me; you have defined your
"lead_shape" function incorrectly, such that no sites are added to your
"dn_lead" system.


Happy Kwanting,

Joe



signature.asc
Description: OpenPGP digital signature