Re: [Kwant] (no subject)

2019-11-04 Thread Abbout Adel 
Dear khani,

I would like to correct my answer and complete  it by answering your
request.

I stated that for an infinite system the conductance is infinite and one
should study something like T(E)/W where W is the width of the system which
we take to infinity.
This is correct. I stated also that we deduce thatintegral  T(E,ky) dky
should also be infinite! That is not correct. Indeed, as long as the
integral is on one Brillouin zone, it will be finite.
(W/pi is the number of Brillouin zones in the finite case!)
I also mentioned that the program I shared ls last time does not wraparound
the scattering region correctly: the reason is that I did a stupid
question. I forgot to put sys[lat.neighbors()]=-1  !
(although the system is one site, we need to put the hoppings since it has
a translational symmetry. )
Below, you can find the correct   program which verifies the analytical
formula I proposed. You can also find how to put a potential in the
scattering region.
The generalization to Graphene is straightforward.

I hope this helps
Adel



from numpy import array,linspace,sqrt,sin,cos,pi,exp,trapz,arccos
import kwant
from matplotlib import pyplot
lat=kwant.lattice.square()
sym1=kwant.TranslationalSymmetry((0,1))
sys=kwant.Builder(sym1)
def sys_pot(site,vg):
return vg
sys[lat(0,0)]=sys_pot
sys[lat.neighbors()]=-1
# in the following put (0,1) first like for sys.
sym2=kwant.TranslationalSymmetry((0,1),(1,0))
lead=kwant.Builder(sym2)
def lead_pot(site,vl):
return vl
lead[lat(0,0)]=lead_pot
lead[lat.neighbors()]=-1

sys=kwant.wraparound.wraparound(sys,coordinate_names='yxz')
lead=kwant.wraparound.wraparound(lead,keep=1,coordinate_names='yxz')
sys.attach_lead(lead)
sys.attach_lead(lead.reversed())

def Trans(E,vg):
transmission=[]
Ky=linspace(0,pi,100)
params={"k_y": 0,"vg":vg, "vl":0}
vg,vl=0,0
for ky in Ky :
params["k_y"]=ky

Sm=kwant.smatrix(sys.finalized(),E,params=params)

transmission.append(Sm.transmission(0,1))
return trapz(transmission,Ky)

Energies=linspace(-2,2,20)
Result=[]   #for vg=0
for E in Energies:
Result.append(Trans(E,0))
#Result.append(Trans(E,0.3)) #if you want vg=0.3

pyplot.plot(Energies,Result,'ro')
# theoretical expected transmission
def T(E):
return arccos(-1+abs(E/2))
Energies=linspace(-2,2,200)
pyplot.plot(Energies,T(Energies))
pyplot.show()
###


On Sun, Oct 27, 2019 at 2:42 PM Abbout Adel  wrote:

> Dear Hosein,
>
> Let me give you my understanding of the problem and then give you an
> example with kwant. I would love to hear a feedback on this matter.
> The best way to deal with this problem is to start with the easiest case:
> 2D square lattice.
> The transmission of an infinite  2D  system is the limit of a waveguide
> with width W when W goes to infinity.
> but we know that the transmission, in this case, is: 2*W/lambda, where
> lambda is the electron wavelength.
>
> So you can see that the conductance (transmission) goes to infinity when
> W---> infinity!
> So, it is better to talk about conductivity rather than conductance
> (conductivity =conductance/W)
>
> I remember that I calculated this for a square lattice during my thesis
> and found:
> T(E)/W =-1/pi   arcos(-1+abs(E/2)) (onsite
> potential v=0)
>
> That means if you do just T(E)=0Intgral T(E, ky) dky , it should diverge
> in 2D!
>
> Let us go back to kwant and see why we get a finite result with a program
> as yours.
> For a uniform 2D system you need to do:
>
> sys=kwant.wraparound.wraparound(sys)
> lead=kwant.wraparound.wraparound(lead, keep=0)
>
> The problem with this is that when you call the scattering matrix:
> Sm=kwant.smatrix(sys.finalized(),energy,[ky])
>
> Only the lead gets the argument ky!  (I don't know why). This means that
> only the lead was wraparounded and got an effective onsite potential (-2t
> cos(ky)) but not the system!
> I checked this analytically for the transmission and found T(E=-1,
> ky=pi/4)=(1+2*sqrt(2))/(3+2*sqrt(2)) which is the same as obtained by kwant
> :0.6568
>
> This confirms my understanding. (I didn't have enough time to think about
> how to correct this)
>
> Another way of accessing the transmission and avoiding what I described
> above is s just to look at the self energy for each ky and each time you
> have an eigenvalue with an imaginary part you have one conducting mode.
>
>
> Your question about the angle: since you have ky, you find the
> corresponding kx at E from the conducting modes and the angle between them
> is what you are searching.
>
> I hope this helps,
> Regards
>
> from numpy import array,linspace,sqrt,sin,cos,pi,exp
> import kwant
> from matplotlib import pyplot
> lat=kwant.lattice.square()
> sym1=kwant.TranslationalSymmetry((0,1))
> sys=kwant.Builder(sym1)
>
> sys[lat(0,0)]=0

Re: [Kwant] (no subject)

2019-10-27 Thread Abbout Adel 
Dear Hosein,

Let me give you my understanding of the problem and then give you an
example with kwant. I would love to hear a feedback on this matter.
The best way to deal with this problem is to start with the easiest case:
2D square lattice.
The transmission of an infinite  2D  system is the limit of a waveguide
with width W when W goes to infinity.
but we know that the transmission, in this case, is: 2*W/lambda, where
lambda is the electron wavelength.

So you can see that the conductance (transmission) goes to infinity when
W---> infinity!
So, it is better to talk about conductivity rather than conductance
(conductivity =conductance/W)

I remember that I calculated this for a square lattice during my thesis and
found:
T(E)/W =-1/pi   arcos(-1+abs(E/2)) (onsite
potential v=0)

That means if you do just T(E)=0Intgral T(E, ky) dky , it should diverge in
2D!

Let us go back to kwant and see why we get a finite result with a program
as yours.
For a uniform 2D system you need to do:

sys=kwant.wraparound.wraparound(sys)
lead=kwant.wraparound.wraparound(lead, keep=0)

The problem with this is that when you call the scattering matrix:
Sm=kwant.smatrix(sys.finalized(),energy,[ky])

Only the lead gets the argument ky!  (I don't know why). This means that
only the lead was wraparounded and got an effective onsite potential (-2t
cos(ky)) but not the system!
I checked this analytically for the transmission and found T(E=-1,
ky=pi/4)=(1+2*sqrt(2))/(3+2*sqrt(2)) which is the same as obtained by kwant
:0.6568

This confirms my understanding. (I didn't have enough time to think about
how to correct this)

Another way of accessing the transmission and avoiding what I described
above is s just to look at the self energy for each ky and each time you
have an eigenvalue with an imaginary part you have one conducting mode.


Your question about the angle: since you have ky, you find the
corresponding kx at E from the conducting modes and the angle between them
is what you are searching.

I hope this helps,
Regards

from numpy import array,linspace,sqrt,sin,cos,pi,exp
import kwant
from matplotlib import pyplot
lat=kwant.lattice.square()
sym1=kwant.TranslationalSymmetry((0,1))
sys=kwant.Builder(sym1)

sys[lat(0,0)]=0

sym2=kwant.TranslationalSymmetry((1,0),(0,1))
lead=kwant.Builder(sym2)
lead[lat(0,0)]=0
lead[lat.neighbors()]=-1

sym3=kwant.TranslationalSymmetry((-1,0),(0,1))
lead2=kwant.Builder(sym3)
lead2[lat(0,0)]=0
lead2[lat.neighbors()]=-1
sys=kwant.wraparound.wraparound(sys)
lead=kwant.wraparound.wraparound(lead,keep=0)
sys.attach_lead(lead)
sys.attach_lead(lead.reversed())



kwant.plot(sys)
def Trans(E):
transmission=[]
Ky=linspace(0,2*pi,100)
for ky in Ky :
Sm=kwant.smatrix(sys.finalized(),energy,[ky])
transmission.append(Sm.transmission(0,1))
return trapz(transmission,Ky)


Energies=linspace(-1.99,-0.0001,30)
Result=[]
for energy in Energies:
Result.append(Trans(energy))

pyplot.plot(Energies,Result)
pyplot.show()




On Fri, Oct 25, 2019 at 12:50 PM Khani Hosein 
wrote:

> Dear all,
> I am using kwant to calculate the transmission of a bulk system. I use "
> sys = kwant.wraparound.wraparound(sys) and  lead =
> kwant.wraparound.wraparound(lead, keep=0)" to change my tight-binding
> ribbon system to a bulk system. The transmission is obtained using:
> kys = np.arange(-0.5, 0.5, 0.001)
> for ky in kys:
> smatrix = kwant.smatrix(sys, energy, [ky])
> transmission.append(smatrix.transmission(1, 0))
>
> I want to obtain the Transmission vs the incident angle of an electron,
> and I also need to do the integration ∫T(ky, E) dky or ∫T(θ,E)cosθdθ. I do
> not undstand the ky in the system with kwant.wraparoun, what is the unit of
> ky, what is k here for the system. Can you give me some advices on
> this?Thanks in advance. The definition of incident angle is shown in the
> attached Fig.1， the band structure for my bulk system is shown in attached
> Fig.2 and the Tranmission vs ky is shown in Fig.3 (I do not know how the
> set the values for ky).
> Regards,
> Hosein Khani
>
>
>

-- 
Abbout Adel

Re: [Kwant] (no subject)

2018-08-10 Thread Abbout Adel 
Dear Kumar,

To do so you need to put a potential in all the leads (you can fix one of
them  to zero) and then use the Conductance matrix.

You will search for V1, V2, (V3=0) such that I1=-I2 and I3=0

(I1,-I1,0)^T  = G   (V1,V2,0)^T

T is the transpose, G is the Conductance matrix

G is singular so to solve this equation, you invert a submatrix of G



def get_Voltage(energy, params, I1=1):

SM=kwant.smatrix(sysf,energy=energy,params=params))
G=SM.conductance_matrix()
M=linalg.inv(G[0:2,0:2])
V=M@array([I1,-I1]).T
 return V


I hope this helps.
Adel

On Fri, Aug 10, 2018 at 9:31 AM, Rohit Kumar Srivastav <
rohit.srivas...@snu.edu.in> wrote:

> I want to apply a voltage on the stadium with three leads in such a way
> two leads constitute the same current and one lead constitute zero current
> (Its means voltage is zero in this lead always) for calculating the
> scattering matrix.
> Please, suggest me what is the code of voltage. (How to apply the voltage)
> Best,
> Rohit
>



-- 
Abbout Adel

Re: [Kwant] (no subject)

2017-03-09 Thread enrique colomes 
Hello,


Definetively I didn't help at all with my previous email.


I want to check that the Haldane model in a ribbon is a gapped system in the 
bulk with edge states. For that purpose, I have modified my jcurrent function 
(attached at the end) and now I obtain correctly the current  through all the 
hopping parameters of the system.


I am trying to plot (not necessarily with Kwant) the total local density 
current in the x and y direction. For that I want to define for each (A/B) atom 
the outgoing current as the vector sum of all the bond currents in each atom 
(each atom has 9 relevant neighbors, except the ones in the edges).


Is there any way of doing that in a simple way with the new current function?

Thanks a lot for your time,

Enrique



def jcurrent(syst):
wfs = kwant.wave_function(syst, energy=0.01)
J = kwant.operator.Current(syst)
wf=wfs(0)[0]
j = J(wf)

#for (i,j), current in zip(syst.graph, j):
   #print('current from site {} to {} is {}'.format(syst.sites[i].tag, 
syst.sites[j].tag, current))




El 07/03/2017 a las 9:25, Joseph Weston escribió:

Hi,



I am using kwant1.3, I have read the new two threads [1] and [2], and I
have read [3]. However, I am still unable to compute and plot the local
density current for my system in the scattering region.


In general you should aim to give more information about precisely what
it is that is not working. Please try and phrase your problem as "I am
trying to do X; I expect the output of this code to be Y, but I receive Z".

Just stating that you are unable to do something does not give anyone
wishing to help enough information to be able to do so.

That being said, I did notice that in the following function:



def jcurrent(sys):

   J=kwant.operator.Current(sys)
   wfs = kwant.wave_function(sys)
   wf= wfs[0]
   j=J(wf)
   kwant.plotter.map(sys, j)


You attempt to plot the output of a kwant operator using
'kwant.plotter.map'. This will raise an error because
'kwant.operator.Current's return the current defined on all *hoppings*
of the system, whereas 'kwant.plotter.map' only knows how to plot
quantities defined on all *sites* in the system. We are currently
working on adding a function 'kwant.plotter.vector_map' that will be
able to plot quantities defined on hoppings (such as the current),
however this is presently not available in Kwant (not even the
development version).

I already responded to Sverre Gulbrandsen overe here [1] where I give
some examples of using current operators. Hopefully that will be enough
to get you started. There will be a full tutorial on operators when
Kwant 1.3 is released.

Happy Kwanting,

Joe


[1]: https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01150.html

Re: [Kwant] (no subject)

2017-03-07 Thread Joseph Weston 
Hi,

> I am using kwant1.3, I have read the new two threads [1] and [2], and I
> have read [3]. However, I am still unable to compute and plot the local
> density current for my system in the scattering region.

In general you should aim to give more information about precisely what
it is that is not working. Please try and phrase your problem as "I am
trying to do X; I expect the output of this code to be Y, but I receive Z".

Just stating that you are unable to do something does not give anyone
wishing to help enough information to be able to do so.

That being said, I did notice that in the following function:

> def jcurrent(sys):
>
>J=kwant.operator.Current(sys)
>wfs = kwant.wave_function(sys)
>wf= wfs[0]
>j=J(wf)
>kwant.plotter.map(sys, j)

You attempt to plot the output of a kwant operator using
'kwant.plotter.map'. This will raise an error because
'kwant.operator.Current's return the current defined on all *hoppings*
of the system, whereas 'kwant.plotter.map' only knows how to plot
quantities defined on all *sites* in the system. We are currently
working on adding a function 'kwant.plotter.vector_map' that will be
able to plot quantities defined on hoppings (such as the current),
however this is presently not available in Kwant (not even the
development version).

I already responded to Sverre Gulbrandsen overe here [1] where I give
some examples of using current operators. Hopefully that will be enough
to get you started. There will be a full tutorial on operators when
Kwant 1.3 is released.

Happy Kwanting,

Joe


[1]: https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01150.html


signature.asc
Description: PGP signature

Re: [Kwant] (no subject)

2017-03-07 Thread Sergey Slizovskiy 

Dear Enrique,
  Below is a bit of code  (basically, it's not mine, it was also taken 
from Kwant discussions), that plots the current density.

Feel free to improve it.
Best wishes,
Sergey

from scipy.sparse import spdiags
def current2d(fsys,lead=0, params=None, h=None, psi=None):
x, y = numpy.array([s.pos for s in fsys.sites]).T
n = len(x)
xd = spdiags(x, 0, n, n)
yd = spdiags(y, 0, n, n)

if h is None:
h = fsys.hamiltonian_submatrix(sparse=True, args=params)
Jx = h * xd - xd * h
Jy = h * yd - yd * h

if psi is None:
wf = kwant.wave_function(fsys, 0, args=params)
#psi = numpy.vstack([wf(l) for l in leads])
psi = wf(lead)
jx = numpy.zeros_like(psi, dtype=float)
jy = numpy.zeros_like(psi, dtype=float)
j2 = numpy.zeros(len(fsys.sites), dtype=float)

for f in xrange(len(psi)):
psi_d_conj = spdiags(psi[f].conj(), 0, n, n)
jx[f] = -2*numpy.imag( psi_d_conj  * (Jx * psi[f]) )
jy[f] = -2*numpy.imag( psi_d_conj  * (Jy * psi[f]) )
 ## plot the wave-function for the first mode in the given lead
   # kwant.plotter.map(fsys, numpy.absolute(psi[f]), 
num_lead_cells=10,a=1,vmin=0,show=False)


jxtot=numpy.sum(jx,axis=0)
jytot=numpy.sum(jy,axis=0)
j2=numpy.sqrt(jxtot**2+jytot**2)###add currents from all the 
modes and then compute the square
kwant.plotter.map(fsys, j2, num_lead_cells=10,a=1,vmin=0, 
vmax=0.1,show=False)


return jxtot , jytot



# If we want to plot the current density a bit of coarse graining helps.
from scipy.interpolate import griddata
from scipy.ndimage.filters import gaussian_filter

def smooth(x, y, jx, jy, L, W, sigma=0, downsample=1, dpu=3):

x_grid = numpy.linspace(-L-1., L, 2*L*dpu+1)
y_grid = numpy.linspace(-1, W+1., W*dpu+1)
grid = numpy.array([_.flatten() for _ in 
numpy.meshgrid(x_grid,y_grid)]).T


jx_grid = griddata((x,y), jx, grid, method='nearest')
jx_grid = jx_grid.reshape(len(y_grid), len(x_grid))

dd = dpu*downsample
#x_gr = x_grid[::dd]
jx_gr = gaussian_filter(jx_grid, sigma*dpu)[::dd,::dd]

jy_grid = griddata((x,y), jy, grid, method='nearest')
jy_grid = jy_grid.reshape(len(y_grid), len(x_grid))
#y_gr = y_grid[::dd]
jy_gr = gaussian_filter(jy_grid, sigma*dpu)[::dd,::dd]  ### takes 
each third element


return jx_gr, jy_gr

# the (coarse-grained) current density looks good in a stream-plot:
def plot_stream(jx, jy, ttl='J', cmap='rainbow', vmax=0.03, Wl=40):
fig, ax = pyplot.subplots(figsize=(11, 7))
jval = numpy.sqrt(jx**2 + jy**2)

sy, sx = jx.shape
#xext =
#yext = (sy-1)
print (sx,sy)
xgr = numpy.linspace(xmin, xmax, sx)
ygr = numpy.linspace(0, Ly, sy)

 #   jval_th = numpy.tanh(jval/vmax)
ax.streamplot(xgr, ygr, jx, jy, density=2, arrowsize=1,
   color = numpy.tanh(jx/vmax), cmap=cmap,
   norm = colors.Normalize(-1,1),
   linewidth=numpy.tanh(jval/vmax)*2)
ax.set_title(ttl)

ax.margins(0.05); ax.set_aspect(1);

# I use bars to indicate the position of the leads around the 
scattering

# region, you need to adjust this to your geometry
##ax.broken_barh([(0, -5)] , (0, xext), edgecolor='grey', 
facecolor='grey')
##ax.broken_barh([(yext, 5)] , (0, xext), edgecolor='grey', 
facecolor='grey')
##ax.broken_barh([(0, Wl)] , (0, -5), edgecolor='grey', 
facecolor='grey')
##ax.broken_barh([(yext, -Wl)] , (0, -5), edgecolor='grey', 
facecolor='grey')
##ax.broken_barh([(0, yext)] , (0, xext), edgecolor='grey', 
facecolor='None')








On 07/03/17 16:20, enrique colomes wrote:

Hello,


I was successful in introducing disorder in my Haldane (second nearest
neighbours) system.


I am using kwant1.3, I have read the new two threads [1] and [2], and I
have read [3]. However, I am still unable to compute and plot the local
density current for my system in the scattering region.


I would really appreciate any help. I attach the code.


Once I have it I will attach it, so people can actually use it and
solve future problems.


Regards,


Enrique


[1]
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html

[2]
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html

[3]

https://kwant-project.org/doc/dev/reference/generated/kwant.operator.Current

from scipy.sparse import spdiags
from math import pi, sqrt, tanh, e
import numpy as np
import kwant
import random
import math
import tinyarray
import scipy.sparse.linalg as sla
from matplotlib import pyplot
# Define the graphene lattice
graphene = kwant.lattice.honeycomb()
a, b = graphene.sublattices

L=4
W=4
t=1
tt= 0.03 * t
phi=1*pi/2


def

Re: [Kwant] (no subject)

2017-03-07 Thread enrique colomes 
Hello,


I was successful in introducing disorder in my Haldane (second nearest
neighbours) system.


I am using kwant1.3, I have read the new two threads [1] and [2], and I
have read [3]. However, I am still unable to compute and plot the local
density current for my system in the scattering region.


I would really appreciate any help. I attach the code.


Once I have it I will attach it, so people can actually use it and
solve future problems.


Regards,


Enrique


[1]
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html

[2]
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html

[3]

https://kwant-project.org/doc/dev/reference/generated/kwant.operator.Current



from scipy.sparse import spdiags
from math import pi, sqrt, tanh, e
import numpy as np
import kwant
import random
import math
import tinyarray
import scipy.sparse.linalg as sla


from matplotlib import pyplot
# Define the graphene lattice
graphene = kwant.lattice.honeycomb()
a, b = graphene.sublattices

L=4
W=4
t=1
tt= 0.03 * t
phi=1*pi/2


def make_system(pot=0.000):

   nnn_hoppinga=tt*e**(1j*phi)
   nnn_hoppingb=tt*e**(1j*phi)

    Define the scattering region. 
   def rec(pos):
x, y = pos
return 0 < x < L and 0 < y < W

   sys = kwant.Builder()

   # w: width and pot: potential maximum of the p-n junction
   def potential(site):
   (x, y) = site.pos
   d = y * 2 + x * 3
   return 0

   sys[graphene.shape(rec, (1, 1))] = potential

   # specify the hoppings of the graphene lattice
   nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
   nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a))
   nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b))
sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] =
-t
sys[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
sys[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb


    Define the leads. 
   # left lead
   sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))

sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)])
sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])

   def lead0_shape(pos):
   x, y = pos
   return (0 < y <  W)

   lead0 = kwant.Builder(sym0)
   lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]]
= -t
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb


   # right lead
   sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0)))

sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)])
sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
   def lead1_shape(pos):
   x, y = pos
   return (0 < y <  W)

   lead1 = kwant.Builder(sym1)
   lead1[graphene.shape(lead1_shape, (0, 0))] = pot
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]]
= -t
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb

   sys.attach_lead(lead1)
   sys.attach_lead(lead0)

   for site in sys.leads[0].interface+sys.leads[1].interface:
   sys[site]=0


   return sys

def jcurrent(sys):
   J=kwant.operator.Current(sys)
   wfs = kwant.wave_function(sys)
   wf= wfs[0]
   j=J(wf)
   kwant.plotter.map(sys, j)


def main():

   sys = make_system()

   # Finalize the system.
   sys = sys.finalized()


   # Plot local current
   jcurrent(sys)


# Call the main function if the script gets executed (as opposed to
imported).
# See 
.
if __name__ == '__main__':
   main()

El 09/02/2017 a las 2:04, Joseph Weston escribió:

Hi again,



b) When I asked about computing it directly from Kwant it was because I read in 
another thread that it was going to be in Kwant1.3, but I was not sure if it 
was already avalaible. I will "fight" then with the code and will try to use 
your guidelines Joseph. Definetively doing it with neighbors will be much 
faster than itering the whole system. Will this new version be ready soon? :-)



Well it is already available in the 'development' version of

Re: [Kwant] (no subject)

2017-03-07 Thread enrique colomes 
Hello,


I was successful in introducing disorder in my Haldane (second nearest
neighbours) system.


I am using kwant1.3, I have read the new two threads [1] and [2], and I
have read [3]. However, I am still unable to compute and plot the local
density current for my system in the scattering region.


I would really appreciate any help. I attach the code.


Once I have it I will attach it, so people can actually use it and
solve future problems.


Regards,


Enrique


[1]
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html

[2]
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html

[3]

https://kwant-project.org/doc/dev/reference/generated/kwant.operator.Current



>from scipy.sparse import spdiags
>from math import pi, sqrt, tanh, e


import numpy as np
import kwant
import random
import math
import tinyarray
import scipy.sparse.linalg as sla


>from matplotlib import pyplot


# Define the graphene lattice
graphene = kwant.lattice.honeycomb()
a, b = graphene.sublattices

L=4
W=4
t=1
tt= 0.03 * t
phi=1*pi/2


def make_system(pot=0.000):

   nnn_hoppinga=tt*e**(1j*phi)
   nnn_hoppingb=tt*e**(1j*phi)

    Define the scattering region. 
   def rec(pos):
x, y = pos
return 0 < x < L and 0 < y < W

   sys = kwant.Builder()

   # w: width and pot: potential maximum of the p-n junction
   def potential(site):
   (x, y) = site.pos
   d = y * 2 + x * 3
   return 0

   sys[graphene.shape(rec, (1, 1))] = potential

   # specify the hoppings of the graphene lattice
   nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
   nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a))
   nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b))
sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] =
-t
sys[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
sys[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb


    Define the leads. 
   # left lead
   sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))

sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)])
sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])

   def lead0_shape(pos):
   x, y = pos
   return (0 < y <  W)

   lead0 = kwant.Builder(sym0)
   lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]]
= -t
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb


   # right lead
   sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0)))

sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)])
sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
   def lead1_shape(pos):
   x, y = pos
   return (0 < y <  W)

   lead1 = kwant.Builder(sym1)
   lead1[graphene.shape(lead1_shape, (0, 0))] = pot
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]]
= -t
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb

   sys.attach_lead(lead1)
   sys.attach_lead(lead0)

   for site in sys.leads[0].interface+sys.leads[1].interface:
   sys[site]=0


   return sys

def jcurrent(sys):
   J=kwant.operator.Current(sys)
   wfs = kwant.wave_function(sys)
   wf= wfs[0]
   j=J(wf)
   kwant.plotter.map(sys, j)


def main():

   sys = make_system()

   # Finalize the system.
   sys = sys.finalized()


   # Plot local current
   jcurrent(sys)


# Call the main function if the script gets executed (as opposed to
imported).
# See 
.
if __name__ == '__main__':
   main()

El 09/02/2017 a las 2:04, Joseph Weston escribió:

Hi again,



b) When I asked about computing it directly from Kwant it was because I read in 
another thread that it was going to be in Kwant1.3, but I was not sure if it 
was already avalaible. I will "fight" then with the code and will try to use 
your guidelines Joseph. Definetively doing it with neighbors will be much 
faster than itering the whole system. Will this new version be ready soon? :-)



Well it is already available in the 'development'

Re: [Kwant] (no subject)

2017-02-08 Thread Abbout Adel 
Dear Enrique,

What I see in your program is that you put a potential
''edge_potential=10'' to help you to color the sites at the edge where you
want to put disorder. In your case, the width is uniform. So you can
achieve your goal  just by putting a constraint on the coordinates of your
sites: Ex:  W-1>site.pos[1]>1.

The most interesting case, is when the shape of your system is arbitrary.
To choose the edge , you do:

for site in sys.expand(graphene.shape(rec, (0, 1))):
if sys.degree(site)<9 :  sys[site]=edge_potential

Please, notice that I used <9 and not <=9 as you did in your program (this
is why everything is green in your plot).

To get rid of the sites at the edge, in front of the leads, just do, after
attaching the leads :

for site in sys.leads[0].interface+sys.leads[1].interface:
   sys[site]=0


This way, your program will work for arbitrary shape, and you can think how
to put disorder after this.
I hope that this helps.
Adel

On Wed, Feb 8, 2017 at 5:31 AM, enrique colomes  wrote:

> Hello,
>
> My name is Enrique and first of all, thanks for the developping of such a
> useful tool.
>
> My question is related to that one of the edge potentials. I am trying to
> reproduce the Haldane model with disorder in a rectangular scattering
> region. My goal is to introduce edge disorder as well as compute the local
> current density.
>
> I have two problems:
>
> a) Regarding the edge disorder, I cannot remove the disorder in the edges
> with the leads (I saw your recent previous discussion)
>
> b) Is it possible already to compute directly from Kwant the local current
> density? If yes, how?
>
> Thanks again,
>
> Enrique
>
>
> from __future__ import division # so that 1/2 == 0.5, and not 0
>
> from scipy.sparse import spdiags
> from math import pi, sqrt, tanh, e
> import numpy as np
> import kwant
> import random
> import math
>
> # For computing eigenvalues
> import scipy.sparse.linalg as sla
>
> # For plotting
> from matplotlib import pyplot
>
> # Define the graphene lattice
> graphene = kwant.lattice.honeycomb()
> a, b = graphene.sublattices
>
>
> edge_potential=10
>
> def make_system(L=10, W=10, pot=0.001):
>
> t=1
> tt= 0.04 * t
> phi=1*pi/2
> edge=1  # 1 if there is edge disorder
> bulk=0  # 1 if there is bulk disorder
> nnn_hoppinga=tt*e**(1j*phi)
> nnn_hoppingb=tt*e**(1j*phi)
>
>  Define the scattering region. 
> def rec(pos):
>  x, y = pos
>  return 0 < x < L and 0 < y < W
>
> sys = kwant.Builder()
>
> # w: width and pot: potential maximum of the p-n junction
> def potential(site):
> (x, y) = site.pos
> d = y * 2 + x * 3
> return 0
>
> sys[graphene.shape(rec, (1, 1))] = potential
>
> # specify the hoppings of the graphene lattice
> nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
> nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a))
> nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b))
> sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t
> sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] 
> = -nnn_hoppinga
> sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] 
> = -nnn_hoppingb
>
>
> # Modify the scattering region
> if edge==1:
> # change the values of the potential at the edge
>  for site in sys.expand(graphene.shape(rec, (0, 1))):
> if sys.degree(site)<=9 and abs(site.pos[0])!=L-4:  
> #abs(site.pos[0])!=19 execludes the interfaces with the leads
> sys[site]=edge_potential
>
> elif bulk==1:
> # change the values of the potential at the edge
>  for site in sys.expand(graphene.shape(rec, (0, 1))):
> if sys.degree(site)==9 and abs(site.pos[0])!=19:  
> #abs(site.pos[0])!=19 execludes the interfaces with the leads
> sys[site]=edge_potential
>
>
>  Define the leads. 
> # left lead
> sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
>
> sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)])
> sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
>
> def lead0_shape(pos):
> x, y = pos
> return (0 < y <  W)
>
> lead0 = kwant.Builder(sym0)
> lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
> lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = 
> -t
> lead0[[kwant.builder.HoppingKind(*hopping) for hopping in 
> nnn_hoppings_a]] = -nnn_hoppinga
> lead0[[kwant.builder.HoppingKind(*hopping) for hopping in 
> nnn_hoppings_b]] = -nnn_hoppingb
>
>
> # right lead
> sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0)))
>
> sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)])
> sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
> def lead1_shape(pos):
> x, y = pos
>

Re: [Kwant] (no subject)

2017-02-08 Thread Joseph Weston 
Hi,

> a) Regarding the edge disorder, I cannot remove the disorder in the edges 
> with the leads (I saw your recent previous discussion)

I am not sure what you mean by this. From what I can see from your model
there is no disorder in your leads:

lead0 = kwant.Builder(sym0)
lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = 
-t 
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] 
= -nnn_hoppinga
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] 
= -nnn_hoppingb

Also, you do not need to manually enumerate all the hopping kinds like this.
Lattices have a method 'neighbors' that returns the hopping kinds for
neighboring sites in the lattice. You could reduce the above code to:

lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
lead0[graphene.neighbors()] = -t
lead0[(kind.family_a == a for kind in graphene.neighbors(2))] = 
-nnn_hoppinga
lead0[(kind.family_a == b for kind in graphene.neighbors(2))] = 
-nnn_hoppingb

which would allow you to remove the lines where you explicitly enumerate the 
hoppings:

nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a))
nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b))


> Is it possible already to compute directly from Kwant the local current 
> density? If yes, how?

This will be available in Kwant version 1.3. I see that at the start of your
script that you have the line 'from __future__ import division' which tells
me that you are still using Python 2. Since version 1.2 Kwant supports Python 3
only, and the 1.1.x versions will only receive bug fixes, no new features.
For this reason I would recommend using Python 3 unless you have a good
reason not to.

Coming back to your question, while calculating currents will be available
in Kwant 1.3, if your model only has 1 orbital per site, it is probably
just as easy to calculate the current manually. The current flowing
from site 'j' to site 'i' is written:

J_ij = 2 Im{(ψ_i)* H_ij ψ_j}

where  'ψ_i' is the component of the wave function on site 'i', '*' is
complex conjugation, and 'H_ij' is the Hamiltonian matrix element
between sites 'i' and 'j'. In Kwant this could be computed like:

psi = kwant.wave_function(syst)(0)[0]  # scattering state 0 from lead 0
current = np.array([2 * psi[i].conjugate() * syst.hamiltonian(i, j) * psi[j]
for i, j in syst.graph])
current = current.imag

In Kwant version 1.3 there will also be functionality for plotting such
quantities defined on hoppings.


Hope that helps,

Joe


signature.asc
Description: PGP signature