Dear everyone,
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.
Thank you for any hint! Best, B
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 = 2 #Length x-direction the length is a bit useless in this code W = 6 #Width y-direction H = 6 #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=25 #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 + heaviside(i,0)*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+V0 #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, 1.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
# Plot transmission plt.plot(energies, Ttot) plt.show()
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
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 christoph.groth@cea.fr wrote:
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
Dear Adel and Krzysiek,
1. thank you for taking time to look into my code and answering my questions in such detail. Sorry also for the late reply. Finally I did not use wraparound. As suggested by Krzysiek, starting from a 3D Hamiltonian in k-space I rewrote the code into a 2D real space model, keeping k_y as a parameter in the leads (and everywhere). The sytem is finite in z direction and x is the transport direction. Then I let kwant obtain the transmission for each k_y and integrated over k_y to obtain the total transmission T_tot(E)=\int dk/(2\pi)T(k_y,E). I hope this is ok.
I am not so sure why I should use wraparound in my situation when I have such a “simple” translational invariance. In the discussion archive I found examples where someone had repeating structures perpendicular to the transport direction, for example like periodical slits/holes in a potential barrier. Probably this is what wraparound is made for. Or maybe also to avoid rebuilding the system for each ky?
2. Just as a remark on a little problem I experienced and would like to share. As I wanted to create bigger slabs I parallelised the loop I used to go through the k_y values. However the code turned out to not benefit at all from this, due to some, as far as I understand, intrinsic parallelisation in the libraries that kwant uses. For any one who is interested, this can be shut down by
os.environ["OMP_NUM_THREADS"] = "1"
then each kwant job runs on one core only and the speed increases linearly with the amount of cores again :).
3. One more quick question please. I would like to plot something like E(ky,kz=0) for my system, and would like to compare it to the exact diagonalisation I did independently of kwant for the uncoupled leads/the bulk on each side of the interface. I found the snippet https://kwant-project.org/doc/1/reference/generated/kwant.physics.Bands#kwan...
bands = kwant.physics.Bands(some_syst) momenta = numpy.linspace(-numpy.pi, numpy.pi, 101) energies = [bands(k) for k in momenta] pyplot.plot(momenta, energies) pyplot.show()
so for my case I can use this piece of code for each ky seperately setting momenta to kz=array([0]) or are there other ways that I overlooked when going through the kwant discussion page/documentation?
Thank you again Adel and dziękuję Krzysiek!
With best regards, Bogusz
On 8 Jun 2019, at 20:51, Abbout Adel abbout.adel@gmail.com wrote:
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:
- 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).
- 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 <christoph.groth@cea.fr mailto:christoph.groth@cea.fr> wrote: 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