Sorry, I just realized I made a silly mistake.

I should replace
sym_lead_vertical = kwant.TranslationalSymmetry((0,-1,0))

by
sym_lead_vertical = kwant.TranslationalSymmetry((0,1,0))

although I don't understand why it is the case.

I apologize for inconvenience caused.

Best,
T.C.

On 20 July 2015 at 00:03, T.C. Wu <[email protected]> wrote:

> Dear all,
>
> I am a beginner in kwant and I am trying to reproduce the quantum Hall
> effect demostrated in
>
> http://nbviewer.ipython.org/github/topocm/topocm_content/blob/master/w3_pump_QHE/Laughlinargument.ipynb
>
> But in my case, I modified the sample code a little bit such that a 3D
> pill-box like system with thickness of 1 atomic layer is built (see
> attached figure). I changed nothing but replaced (xs,ys) by (xs,ys,zs) and
> so on. However, the Hall conductance simulated in my system is completely
> different from the 2D case ( no plateaus formed, see the attached figure).
>
>
> Could you please give me some advice?
>
>
> # My code#
>
> import numpy as np
> import matplotlib
> import kwant
> from matplotlib import pyplot as plt
>
> L = 100
> W = 60
> H = 1
> w_lead = W-2
> w_vert_lead = L/2-1
> t = 1.0
> mu = 0.3
> mu_lead = 0.3
> Bs = np.linspace(.02, .15, 200)
>
> class SimpleNamespace(object):
>     def __init__(self, **kwargs):
>         self.__dict__.update(kwargs)
>
>
> ## Quantum hall bar codes
> def qhe_hall_bar(L=50, W=10,H=2, w_lead=10, w_vert_lead=None):
>     """Create a hall bar system.
>
>     Square lattice, one orbital per site.
>     Returns finalized kwant system.
>
>     Arguments required in onsite/hoppings:
>         t, mu, mu_lead
>     """
>
>     L = 2*(L//2)
>     W = 2*(W//2)
>     H = 2*(H//2)
>     w_lead = 2*(w_lead//2)
>     if w_vert_lead is None:
>         w_vert_lead = w_lead
>     else:
>         w_vert_lead = 2*(w_vert_lead//2)
>
>     # bar shape
>     def bar(pos):
>         (x, y,z) = pos
>         return (x >= -L/2 and x <= L/2) and (y >= -W/2 and y <= W/2) and
> (z>= -H/2 and z <= H/2)
>
>     # Onsite and hoppings
>     def onsite(site, par):
>         (x,y,z) = site.pos
>         return 4*par.t - par.mu
>
>     def hopping_Ax(target, source, par):
>         xt, yt, zt = target.pos
>         xs, ys, zs = source.pos
>         return -par.t * np.exp(-0.5j * par.B * (xt + xs) * (yt - ys))
>
>     def make_lead_hop_y(x0):
>         def hopping_Ay(target, source, par):
>             print (x0)
>             xt, yt, zt = target.pos
>             xs, ys, zs = source.pos
>             return -par.t * np.exp(-1j * par.B * x0 * (yt - ys))
>         return hopping_Ay
>
>     def lead_hop_vert(target, source, par):
>         return -par.t
>
>     # Building system
>     lat = kwant.lattice.general(np.identity(3))
>     sys = kwant.Builder()
>     print('sys built')
>
>     sys[lat.shape(bar, (0, 0,0))] = onsite
>     sys[lat.neighbors()] = hopping_Ax
>     print('sys hopping done')
>     #kwant.plot(sys)
>     # Attaching leads
>     sym_lead = kwant.TranslationalSymmetry((-1, 0,0))
>
>
>     def lead_shape(pos):
>         (x, y,z) = pos
>         print('lead_shape')
>         return (-w_lead / 2 <= y <= w_lead / 2) and (z>= -H/2 and z <= H/2)
>
>     lead_onsite = lambda site, par: 4 * par.t - par.mu_lead
>
>     sym_lead_vertical = kwant.TranslationalSymmetry((0,-1,0))
>     lead_vertical1 = kwant.Builder(sym_lead_vertical)
>     lead_vertical2 = kwant.Builder(sym_lead_vertical)
>     print('vert lead done')
>     def lead_shape_vertical1(pos):
>         (x, y,z) = pos
>         return (-L/4 - w_vert_lead/2 <= x <= -L/4 + w_vert_lead/2) and
> (H/2 >= z>= -H/2)
>
>     def lead_shape_vertical2(pos):
>         (x, y,z) = pos
>         return (+L/4 - w_vert_lead/2 <= x <= +L/4 + w_vert_lead/2) and
> (H/2 >= z>= -H/2)
>
>     lead_vertical1[lat.shape(lead_shape_vertical1, (-L/4 , 0,0))] =
> lead_onsite
>     print('onsite')
>     lead_vertical1[lat.neighbors()] = lead_hop_vert
>     lead_vertical2[lat.shape(lead_shape_vertical2, (L/4,0, 0))] =
> lead_onsite
>     lead_vertical2[lat.neighbors()] = lead_hop_vert
>
>     print('lead defined')
>
>     sys.attach_lead(lead_vertical1)
>     sys.attach_lead(lead_vertical2)
>
>     sys.attach_lead(lead_vertical1.reversed())
>     sys.attach_lead(lead_vertical2.reversed())
>     lead = kwant.Builder(sym_lead)
>     lead[lat.shape(lead_shape, (-1, 0,0))] = lead_onsite
>     lead[lat.neighbors()] = make_lead_hop_y(-L/2)
>
>     sys.attach_lead(lead)
>
>     lead = kwant.Builder(sym_lead)
>     lead[lat.shape(lead_shape, (-1, 0,0))] = lead_onsite
>     lead[lat.neighbors()] = make_lead_hop_y(L/2)
>
>     sys.attach_lead(lead.reversed())
>     print('lead attached')
>     return sys.finalized()
>
> def calculate_sigmas(G):
> #     smat = smat.data
> #     smat = abs(smat)**2
>     # reduce by one dimension G -> G[temp, temp]
>     temp = [0, 1, 3, 4, 5]
>     G = G[temp, :]
>     G = G[:, temp]
>     # invert R = G^-1
>     # find out whether it is a numpy object
>     r = np.linalg.inv(G)
>     # Voltages follow: V = R I[temp]
>     V = r.dot(np.array([0, 0, 0, -1, 1]))
>     # Completely solved the six terminal system.
>     # Consider the 2x2 conductance now: Use I = sigma U
>     E_x = V[1] - V[0]
>     E_y = V[1] - V[3]
>     #formula above
>     sigma_xx = E_x / (E_x**2 + E_y**2)
>     sigma_xy = E_y / (E_x**2 + E_y**2)
>     return sigma_xx, sigma_xy ,-1*E_x,-1*E_y
> print('hi')
> par = SimpleNamespace(t=t, mu=mu, mu_lead=mu_lead, B=None)
> sys = qhe_hall_bar(L=L, W=W,H=H ,w_lead=w_lead, w_vert_lead=w_vert_lead)
>
> kwant.plot(sys)
> sigmasxx = []
> sigmasxy = []
> rxx = []
> rxy = []
> #Bs = np.linspace(.02, .15, 200)
>
> for par.B in Bs:
>     s = kwant.smatrix(sys, args=[par])
>     G = np.array([[s.transmission(i, j) for i in range(len(sys.leads))]
>                   for j in range(len(sys.leads))])
>     G -= np.diag(np.sum(G, 0))
>     sigma = calculate_sigmas(G)
>     sigmasxx.append(sigma[0])
>     sigmasxy.append(sigma[1])
>     rxx.append(sigma[2])
>     rxy.append(sigma[3])
>
> fig = plt.figure()
> ax = plt.subplot(1, 1, 1)
>
> ax.plot(1/Bs, sigmasxx, 'k')
> ax.plot(1/Bs, sigmasxy, 'r');
> ax.set_xticks([])
> ax.set_xlabel('$B^{-1} [a.u.]$')
> ax.set_yticks(range(4))
> ax.set_yticklabels(['${0}$'.format(i) for i in range(4)])
> ax.set_ylabel('$\sigma_{xx}, \sigma_{xy}\,[e^2/h]$')
> plt.show()
>
> plt.plot(Bs,rxx)
> plt.plot(Bs,rxy)
> plt.show()
> #ax.set_ylim([-.1, 3.3]);
>
>

Reply via email to