Dear Bogusz,
I checked your program and it seems fine and the result should be 'almost'
correct.
However, I would like to make some comments and corrections:
1) You are using wraparound and your system if uniform in the y-direction.
No need to use a width W. A one unit cell length is enough.
This may seem at a first glance with no consequences, but it does change
the definition of your ky such that two people using two different widths
will have two different values of T(ky, E) although the integrated value
T(E) should be the same (if integrated on the correct interval in each
case).
2) Let us first be interested on the sum of T(ky, E) over some chosen
continuous modes ky and discretized open modes in z-direction. Let us also
put the same two leads at the same potential and put the L=1. This will
help us in doing analytic calculation.
The energy of each mode is:
Em=6*t-2*t*cos(ky)-2*t*cos(m*pi/(H+1))
There is no inter-mode scattering in your system and the transmission can
be obtained using the Fisher-Lee formula as for a scatterer in 1D:
T(ky,E)=sum_m 1/(1.+(V0/Gamma)**2) with Gamma=sqrt(4-(E-Em)**2).
This formula is tested in the below code.
If you change W>1, it will not work for the reason I mentioned before: some
re normalization of the ky (ky/W) and changes should be done probably. I
tried to figure out how to do it but I failed.
If you put V0 in one lead (as in your program) a quick calculation can give
you the analytical result in the same way as here (you will need Gamma1,
Gamma2)
I hope this helps
Adel
############################
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from numpy import *
import kwant
import sys
#calculate the conductance for a sheet (Dim=2) or a Slab (Dim=3) that is
translational invariant in y direction
a = 1.0 #lattice constant
V0 = 0.2 #potential step size
L = 1 #Length x-direction the length is a bit useless in
this code
W = 1 #Width y-direction
H = 4 #Hight z-direction
t = 1 #hopping parameter
Dim=3 #dimension of the system
# for H=1, Dim=2
# for H>1, Dim=3
#params for integration and plotting
Nky=100 #points in ky space
NE=30 #points in energy
lat = kwant.lattice.cubic(a)
# Define the scattering region
# Infinite potential plane/Slab in y direction
sys = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((0, W, 0))))
for i in range(L):
for j in range(W):
for k in range(H):
# On-site Hamiltonian
sys[lat(i, j, k)] = 2*Dim * t + V0 #potential
step located at the center between the leads
sys[lat.neighbors(1)] = -t
sys = kwant.wraparound.wraparound(sys)
#Left lead
leadL = kwant.Builder(kwant.TranslationalSymmetry((-a, 0, 0), lat.vec((0,
W, 0)))) #trans inv in transport dir and ky
leadL[(lat(0,j,k) for j in range(W) for k in range(H))] = 2*Dim*t
#onsite terms in the left lead
leadL[lat.neighbors(1)] = -t
#Right lead
leadR = kwant.Builder(kwant.TranslationalSymmetry(( a, 0, 0), lat.vec((0,
W, 0))))
leadR[(lat(0,j,k) for j in range(W) for k in range(H))] = 2*Dim*t
#onsite terms in the right lead
leadR[lat.neighbors(1)] = -t
leadL = kwant.wraparound.wraparound(leadL, keep=0) #keeping the
translational invariance in x direction
leadR = kwant.wraparound.wraparound(leadR, keep=0)
sys.attach_lead(leadL)
sys.attach_lead(leadR)
kwant.plot(sys)
sys = sys.finalized()
# Calculation of the total transmission for a given Energy
ky_array = np.linspace(-pi,pi,Nky)
dky=ky_array[1]-ky_array[0]
energies = np.linspace(0.0, 5.0,NE)
Ttot=[]
for energy in energies:
Tky_array=zeros(Nky)
for i in arange(len(ky_array)):
ky=ky_array[i]
smatrix = kwant.smatrix(sys, energy, [ky])
Tky_array[i]=smatrix.transmission(1, 0)
#Ttot.append((2*pi)/(W)*(sum(Tky_array)-0.5*Tky_array[0]-0.5*Tky_array[Nky-1])*dky)
#cheap trapezodial rule to integrate over ky
Ttot.append(sum(Tky_array)) #cheap trapezodial rule to integrate over ky
# Plot transmission
plt.plot(energies, Ttot)
plt.show()
def T(E,qy,W,H,V0,t=1):
T0=0
for m in range(1,H+1):
Em=6*t-2*t*cos(qy/W)-2*t*cos(m*pi/(H+1))
if abs(E-Em)<=2: # condition for open modes
Gamma=sqrt(4-(E-Em)**2)
T0+=1/(1.+(V0/Gamma)**2) ##use Fisher Lee Formula T=Gamma G
Gamma G^\dagger
return T0
Trans=[sum(np.array([T(E,qy,W,H,V0,t=1) for qy in ky_array ] )) for E in
energies]
plt.plot(energies,Trans,'ro')
plt.plot(energies, Ttot)
plt.show()
On Thu, Jun 6, 2019 at 7:03 PM Christoph Groth
Bogusz Bujnowski wrote:
I would like to solve the 3D Schrödinger equation with a 1D potential step. With the material in the mailing list as well as from the tutorial the 2D case works fine.
Extending to 3D gives the same when just taking one lattice site in the z direction (H>1). Increasing the hight of the slab (H>1) changes the minimal energy for which transmission is possible thus my code is wrong.
Did I define and attach the leads correctly with the wraparound in 3D? The short code is below.
From a quick look at your script, the way you are using wraparound seems fine. Unfortunately, I don’t have the time to debug your problem, but here’s a suggestion of what you can do:
You are using a 3D model with two translational symmetries. One of them you convert to a momentum parameter using wraparound, and the other one you use for the leads.
You could simplify your script to use a 2D model, but still use two translational symmetries in the same way as now. This way, you can more easily examine your system by plotting it, and the calculations are quicker.
When plotting your wrapped-around system in 2D, you should see hoppings that "wrap around" (hence the name). This should be equivalent to creating such a Builder by hand (and it should be quite easy to do in your case). So you should be able to debug the problem.
Please let us know when you find a solution. We are particularly interested if you find some problem with Kwant (or its documentation).
Owocnego Kwantowania!
Krzysiek
-- Abbout Adel