From physical reasoning, I would expect that the current through each link in the direction of the current is the same. So I tried to apply the formula (11) from the NJP-introduction to KWANT. But I get a very different looking answers, comparing to scattering matrix approach, as well as current throught different links in the direction of
Thanks a lot for you package which is rather intuitive in using! Could you please help with understanding of the wave function which KWANT gives? I could not understand what kwant.wave_function gives as an output: -- how to distinguish between the modes which are propagating in the leads and which are decaying in the leads? kwant.wave_function gives an answer at any energy. (1)For example, how do get the density in some transport setup? Is this example correct? def Fermi(energy, T, mu): f = 1./(1. + math.exp ( (energy - mu) / T) ) return f def density(sys,energy,lead_nr): wf = kwant.wave_function(sys,energy) return (abs(wf(lead_nr))**2).sum(axis=0) numleads = 2 Tmu = numpy.array([[Tleft,muleft],[Tright,muright]]) #assuming that this range of energies encompasses well the chemical potential range de = 0.05 nenergies = 100 energies=[de * (i - nenergies/2) for i in xrange(nenergies)] for energy in energies: for num_lead in xrange(numleads): d = d + Fermi(energy, Tmu[num_lead]) * density(sys, energy, num_lead) d = d * de (2) Another example which I do not understand properly is about transport through a translationally invariant 2d bar. the current is different. Why are these two answers so different and one is not oscillating, and another is so much depends on energy? (In 1d the comparission works well) def Fermi(energy, T, mu): f = 1./(1. + math.exp ( (energy - mu) / T) ) return f #potentials in the leads muleft = -2. muright = 2. mu = 0 Tleft = 0.1 Tright = 0.1 #potentials and temepratures of all leads muTAll = [] muTAll.append([muleft,Tleft]) muTAll.append([muright,Tright]) energies = [] data = [] nenergies = 300 de = 0.032424 CurrentWave = [] CurrentLandauer = [] t = 1. W = 6 L = 10 numleads = 2 numsite = (L-1)* W sys = BuildSystem(t, W, L,mu,mu,mu) #beginning and end of stm link in terms of program; x=[] for i1 in xrange(numsite): x.append(sys.sites[i1].pos) beg=[] for i1 in xrange(numsite): if x[i1] == [1,1]: beg.append(i1) end=[] for i1 in xrange(numsite): if x[i1] == [2,1]: end.append(i1) print beg, end for imu in xrange(1): sysFinal = BuildSystem(t, W, L, mu, mu, mu) for ienergy in xrange(nenergies): energy = de * ienergy - nenergies * de / 2. energies.append(energy) wf = kwant.wave_function(sysFinal,energy) smatrix = kwant.smatrix(sysFinal, energy) #wave function based calculation which can be used for any correlation function #iterating over number of leads wfInside = numpy.zeros(numsite) #wave function inside of the system for p in xrange( numleads): wave = wf(p) lenwave = wave.shape [0] for ipsi in xrange(lenwave): wfInside = wfInside + Fermi(energy, muTAll[p][1],muTAll[p][0] ) * wave[ipsi] ProbeCurrent = wfInside[beg] * wfInside[end].conjugate() - wfInside[end] * wfInside[beg].conjugate() CurrentWave.append( ProbeCurrent[0].imag) #Getting wave function is done #transmission based calculation t10 = smatrix.transmission(1, 0) data.append(t10) t2 = t10 * (Fermi(energy, muTAll[1][1],muTAll[1][0] ) - Fermi(energy, muTAll[0][1],muTAll[0][0] )) CurrentLandauer.append(t2) plt.figure() plt.plot(energies, CurrentLandauer) plt.plot(energies, CurrentWave) plt.xlabel("energy [t]") plt.ylabel("conductance [e^2/h]") plt.show() def BuildSystem(t, W, L, mu, muleft, muright): sys = kwant.Builder() #DEFINING THE SYSTEM a=1. lat = kwant.lattice.square(a) for i in xrange(L-1): for j in xrange(W): sys[lat(i,j)] = mu if j > 0: sys[lat(i,j), lat(i,j-1)] = -t if i > 0: sys[lat(i,j), lat(i-1,j)] = -t #ATTACHING LEADS #lead is pointing away from scattering region sym_left_lead = kwant.TranslationalSymmetry((-a,0)) left_lead = kwant.Builder(sym_left_lead) for j in xrange(W): left_lead[lat(0,j)] = muleft if j > 0: #modeling large band-width in the leads, that is why 10. t left_lead[lat(0,j), lat(0,j-1)] = - t left_lead[lat(1,j), lat(0,j)] = - t sys.attach_lead(left_lead) sym_right_lead = kwant.TranslationalSymmetry((a,0)) right_lead = kwant.Builder(sym_right_lead) for j in xrange(W): numx = L-2 right_lead[lat(numx,j)] = muright if j > 0: #modeling large band-width in the leads, that is why 10. t right_lead[lat(numx,j), lat(numx,j-1)] = - t right_lead[lat(numx,j), lat(numx+1,j)] = - t sys.attach_lead(right_lead) sysFinal = sys.finalized() t2 = time.clock() return sysFinal
From physical reasoning, I would expect that the current through each link in the direction of the current is the same. So I tried to apply the formula (11) from the NJP-introduction to KWANT. But I get a very different looking answers, comparing to scattering matrix approach, as well as current throught different links in the direction of
Thanks a lot for you package which is rather intuitive in using! Could you please help with understanding of the wave function which KWANT gives? I could not understand what kwant.wave_function gives as an output: -- how to distinguish between the modes which are propagating in the leads and which are decaying in the leads? kwant.wave_function gives an answer at any energy. (1)For example, how do get the density in some transport setup? Is this example correct? def Fermi(energy, T, mu): f = 1./(1. + math.exp ( (energy - mu) / T) ) return f def density(sys,energy,lead_nr): wf = kwant.wave_function(sys,energy) return (abs(wf(lead_nr))**2).sum(axis=0) numleads = 2 Tmu = numpy.array([[Tleft,muleft],[Tright,muright]]) #assuming that this range of energies encompasses well the chemical potential range de = 0.05 nenergies = 100 energies=[de * (i - nenergies/2) for i in xrange(nenergies)] for energy in energies: for num_lead in xrange(numleads): d = d + Fermi(energy, Tmu[num_lead]) * density(sys, energy, num_lead) d = d * de (2) Another example which I do not understand properly is about transport through a translationally invariant 2d bar. the current is different. Why are these two answers so different and one is not oscillating, and another is so much depends on energy? (In 1d the comparission works well) def Fermi(energy, T, mu): f = 1./(1. + math.exp ( (energy - mu) / T) ) return f #potentials in the leads muleft = -2. muright = 2. mu = 0 Tleft = 0.1 Tright = 0.1 #potentials and temepratures of all leads muTAll = [] muTAll.append([muleft,Tleft]) muTAll.append([muright,Tright]) energies = [] data = [] nenergies = 300 de = 0.032424 CurrentWave = [] CurrentLandauer = [] t = 1. W = 6 L = 10 numleads = 2 numsite = (L-1)* W sys = BuildSystem(t, W, L,mu,mu,mu) #beginning and end of stm link in terms of program; x=[] for i1 in xrange(numsite): x.append(sys.sites[i1].pos) beg=[] for i1 in xrange(numsite): if x[i1] == [1,1]: beg.append(i1) end=[] for i1 in xrange(numsite): if x[i1] == [2,1]: end.append(i1) print beg, end for imu in xrange(1): sysFinal = BuildSystem(t, W, L, mu, mu, mu) for ienergy in xrange(nenergies): energy = de * ienergy - nenergies * de / 2. energies.append(energy) wf = kwant.wave_function(sysFinal,energy) smatrix = kwant.smatrix(sysFinal, energy) #wave function based calculation which can be used for any correlation function #iterating over number of leads wfInside = numpy.zeros(numsite) #wave function inside of the system for p in xrange( numleads): wave = wf(p) lenwave = wave.shape [0] for ipsi in xrange(lenwave): wfInside = wfInside + Fermi(energy, muTAll[p][1],muTAll[p][0] ) * wave[ipsi] ProbeCurrent = wfInside[beg] * wfInside[end].conjugate() - wfInside[end] * wfInside[beg].conjugate() CurrentWave.append( ProbeCurrent[0].imag) #Getting wave function is done #transmission based calculation t10 = smatrix.transmission(1, 0) data.append(t10) t2 = t10 * (Fermi(energy, muTAll[1][1],muTAll[1][0] ) - Fermi(energy, muTAll[0][1],muTAll[0][0] )) CurrentLandauer.append(t2) plt.figure() plt.plot(energies, CurrentLandauer) plt.plot(energies, CurrentWave) plt.xlabel("energy [t]") plt.ylabel("conductance [e^2/h]") plt.show() def BuildSystem(t, W, L, mu, muleft, muright): sys = kwant.Builder() #DEFINING THE SYSTEM a=1. lat = kwant.lattice.square(a) for i in xrange(L-1): for j in xrange(W): sys[lat(i,j)] = mu if j > 0: sys[lat(i,j), lat(i,j-1)] = -t if i > 0: sys[lat(i,j), lat(i-1,j)] = -t #ATTACHING LEADS #lead is pointing away from scattering region sym_left_lead = kwant.TranslationalSymmetry((-a,0)) left_lead = kwant.Builder(sym_left_lead) for j in xrange(W): left_lead[lat(0,j)] = muleft if j > 0: #modeling large band-width in the leads, that is why 10. t left_lead[lat(0,j), lat(0,j-1)] = - t left_lead[lat(1,j), lat(0,j)] = - t sys.attach_lead(left_lead) sym_right_lead = kwant.TranslationalSymmetry((a,0)) right_lead = kwant.Builder(sym_right_lead) for j in xrange(W): numx = L-2 right_lead[lat(numx,j)] = muright if j > 0: #modeling large band-width in the leads, that is why 10. t right_lead[lat(numx,j), lat(numx,j-1)] = - t right_lead[lat(numx,j), lat(numx+1,j)] = - t sys.attach_lead(right_lead) sysFinal = sys.finalized() t2 = time.clock() return sysFinal
Dear Mariya, Kwant’s wave_function routine will give you the wave functions inside the scattering region due to all the modes of a lead that are incoming and propagating. Does this answer your question? Christoph
Dear Christoph, Thanks for feedback, but I am still confused. Let me make an example: there are two leads 0 and 1, the lead 0 is filled more than lead 1; from the leads 0 there is an outgoing mode exp(i k x), which is normalized by flux in the leads; the wave function inside of the scattering region is psi = t exp (i k x) + r exp (-i k x) What does the program give if the wave function is called with lead number 0? psi? (and when the leads have couple of modes, then psi[i] are sorted is some comfortable for the program way) or r exp (-i k x)? Best, Mariya
Hi, I'd like to add to Christoph's answer
Thanks for feedback, but I am still confused. Let me make an example: there are two leads 0 and 1, the lead 0 is filled more than lead 1; from the leads 0 there is an outgoing mode exp(i k x), which is normalized by flux in the leads; the wave function inside of the scattering region is psi = t exp (i k x) + r exp (-i k x)
What does the program give if the wave function is called with lead number 0? psi? (and when the leads have couple of modes, then psi[i] are sorted is some comfortable for the program way) or r exp (-i k x)?
Calling the wavefunction with lead number 0 would give you `psi`, the *full scattering state*, in your example. As to the ordering, the kwant documentation says (http://bit.ly/1EziBYZ http://bit.ly/1EziBYZ): The modes appear in the same order as incoming modes in kwant.physics.modes http://kwant-project.org/doc/1.0/reference/generated/kwant.physics.modes.htm... Hopefully this clarifies things. In your previous email you asked about calculating the electronic density in a transport setup and gave an example. From what I can see your example is correct in the case that there are no true bound states in the system. Any bound states will contribute to the density and are not calculated by kwant.wave_function (but they should not contribute to the DC transport). In the second example you ask about the equivalence between calculating a current in the leads using scattering matrices and current in the central part of the system by using the scattering wavefunctions. You are, of course, correct in your statement that the two should give the same result (at least in your example of a bar attached to 2 leads). From an initial glance at your code, it seems that there are at least 2 problems with your wavefunction calculation: 1) You add together the wavefunctions corresponding to different modes, and then calculate the current given by this coherent superposition. What you should actually do is calculate the current due to *each* of the wavefunctions and then add these currents. 2) When you calculate the current you only do so between sites (1, 1) and (2, 1), despite the fact that your system is 6 sites wide. To get the total current flowing through the bar you should, of course, calculate the current between sites (1, i) -> (2, i) for i in (0, 6] and add them all together. Best, Joe
Dear Joseph,
Thanks a lot for clarification.
Best,
Mariya
On 28 April 2015 at 14:55, Joseph Weston
Hi,
I'd like to add to Christoph's answer
Thanks for feedback, but I am still confused. Let me make an example: there are two leads 0 and 1, the lead 0 is filled more than lead 1; from the leads 0 there is an outgoing mode exp(i k x), which is normalized by flux in the leads; the wave function inside of the scattering region is psi = t exp (i k x) + r exp (-i k x)
What does the program give if the wave function is called with lead number 0? psi? (and when the leads have couple of modes, then psi[i] are sorted is some comfortable for the program way) or r exp (-i k x)?
Calling the wavefunction with lead number 0 would give you `psi`, the *full scattering state*, in your example. As to the ordering, the kwant documentation says (http://bit.ly/1EziBYZ):
The modes appear in the same order as incoming modes in kwant.physics.modes http://kwant-project.org/doc/1.0/reference/generated/kwant.physics.modes.htm...
Hopefully this clarifies things.
In your previous email you asked about calculating the electronic density in a transport setup and gave an example. From what I can see your example is correct in the case that there are no true bound states in the system. Any bound states will contribute to the density and are not calculated by kwant.wave_function (but they should not contribute to the DC transport).
In the second example you ask about the equivalence between calculating a current in the leads using scattering matrices and current in the central part of the system by using the scattering wavefunctions. You are, of course, correct in your statement that the two should give the same result (at least in your example of a bar attached to 2 leads). From an initial glance at your code, it seems that there are at least 2 problems with your wavefunction calculation:
1) You add together the wavefunctions corresponding to different modes, and then calculate the current given by this coherent superposition. What you should actually do is calculate the current due to *each* of the wavefunctions and then add these currents. 2) When you calculate the current you only do so between sites (1, 1) and (2, 1), despite the fact that your system is 6 sites wide. To get the total current flowing through the bar you should, of course, calculate the current between sites (1, i) -> (2, i) for i in (0, 6] and add them all together.
Best,
Joe
participants (5)
-
Christoph Groth
-
Joseph Weston
-
Mariya
-
Mariya Medvedyeva
-
Mariya Medvedyeva