Good day, I'm trying to reconstruct something discussed in http://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00432.html. In this answer it is described how to take the wavefunction from a 3D system and plot it in a 2D slice. This turned out to be something different than what the OP was asking for, but that is not so important. It is also mentioned that this has been asked before, but I couldn't find the right search terms. If my question has already been answered, I'd be very grateful for a referral to the relevant discussion! Now, as it turns out, I can't seem to get the above to work. I take some small cylinder, I calculate the wavefunction, but then I get stuck at how to extract for example the y = 0 plane and plot this in 2D. Below I'll post how far I got, perhaps the fix would be rather simple. def onsite(site, p): return 6*p.t def hopping(site0, site1, p): return -p.t def make_system_cyl(a=1, s=10, L = 15): lat = kwant.lattice.general([(0, 0, a), (a, 0, 0), (0, a, 0)]) sys = kwant.Builder() def cyl_shape(pos): x, y, z = pos rsq = y ** 2 + z ** 2 return rsq < s^2 and 0 <= x < L sys[lat.shape(cyl_shape, (1,1,1))] = onsite sys[lat.neighbors()] = hopping def lead_shape(pos): x, y, z = pos rsq = y ** 2 + z ** 2 return rsq < s^2 sym_lead = kwant.TranslationalSymmetry((-a,0,0)) lead = kwant.Builder(sym_lead) lead[lat.shape(lead_shape, (-a,1,1))] = onsite lead[lat.neighbors()] = hopping sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys and I extract the wavefunction at some energy sys = make_system_cyl(a=1, s=10, L = 15) kwant.plot(sys) sysf = sys.finalized() params = SimpleNamespace(t=1) wf = kwant.wave_function(sysf, 1.5, args=[params]) And then I try to extract the relevant parts wavefunction = wf(0) projected_wavefunction = np.empty((wavefunction.size, 1), dtype=np.complex) def in_sheet(index, site): return site.pos[1] == 0 # site in x-z plane sheet_idx_sites = filter(in_sheet, enumerate(sys.sites) for i, (idx, site) in enumerate(sheet_idx_sites): projected_wavefunction[i] = wavefunction[idx] Now this last part I took from the post referenced at the start, and it gives a syntax error. I'm not really sure how to fix it as it looks fine to me. What's wrong there? And more importantly, once this is fixed, how can one then make a kwant.plotter.map type of plot at this x-z plane?
Sorry, the above code has two typo's at the end. It should be sheet_idx_sites = filter(in_sheet, enumerate(sysf.sites)) for i, (idx, site) in enumerate(sheet_idx_sites): projected_wavefunction[i] = wavefunction[idx] This now gives the error that in_sheet() missing 1 required positional argument: 'site' 2017-04-25 16:43 GMT+02:00 John Eaves <johneaves02@gmail.com>:
Good day,
I'm trying to reconstruct something discussed in http://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00432.html. In this answer it is described how to take the wavefunction from a 3D system and plot it in a 2D slice. This turned out to be something different than what the OP was asking for, but that is not so important. It is also mentioned that this has been asked before, but I couldn't find the right search terms. If my question has already been answered, I'd be very grateful for a referral to the relevant discussion!
Now, as it turns out, I can't seem to get the above to work. I take some small cylinder, I calculate the wavefunction, but then I get stuck at how to extract for example the y = 0 plane and plot this in 2D. Below I'll post how far I got, perhaps the fix would be rather simple.
def onsite(site, p): return 6*p.t
def hopping(site0, site1, p): return -p.t
def make_system_cyl(a=1, s=10, L = 15): lat = kwant.lattice.general([(0, 0, a), (a, 0, 0), (0, a, 0)])
sys = kwant.Builder()
def cyl_shape(pos): x, y, z = pos rsq = y ** 2 + z ** 2 return rsq < s^2 and 0 <= x < L
sys[lat.shape(cyl_shape, (1,1,1))] = onsite sys[lat.neighbors()] = hopping
def lead_shape(pos): x, y, z = pos rsq = y ** 2 + z ** 2 return rsq < s^2
sym_lead = kwant.TranslationalSymmetry((-a,0,0)) lead = kwant.Builder(sym_lead)
lead[lat.shape(lead_shape, (-a,1,1))] = onsite lead[lat.neighbors()] = hopping
sys.attach_lead(lead) sys.attach_lead(lead.reversed())
return sys
and I extract the wavefunction at some energy
sys = make_system_cyl(a=1, s=10, L = 15) kwant.plot(sys) sysf = sys.finalized() params = SimpleNamespace(t=1) wf = kwant.wave_function(sysf, 1.5, args=[params])
And then I try to extract the relevant parts
wavefunction = wf(0) projected_wavefunction = np.empty((wavefunction.size, 1), dtype=np.complex)
def in_sheet(index, site): return site.pos[1] == 0 # site in x-z plane
sheet_idx_sites = filter(in_sheet, enumerate(sys.sites)
for i, (idx, site) in enumerate(sheet_idx_sites): projected_wavefunction[i] = wavefunction[idx]
Now this last part I took from the post referenced at the start, and it gives a syntax error. I'm not really sure how to fix it as it looks fine to me. What's wrong there?
And more importantly, once this is fixed, how can one then make a kwant.plotter.map type of plot at this x-z plane?
Hi John, If I understand correctly you are trying to visualize the square magnitude of some wavefunction of a 2D section of a 3D system. This problem has 2 parts: first the relevant wavefunctions components must be extracted from the wavefunction defined over the full (3D) system, and then these components must be plotted, e.g. by plotting a 2D colormap. The first part is relatively straighforward (and indeed you found a solution on the mailing list), however the second part has not been explicitly discussed on the mailing list and is indeed more challenging. First to tackle the mistakes that I saw in your script. In your system creation function you have the line return rsq < s^2 when what you probably want is return rsq < s**2 The caret (^) operator in Python means bitwise XOR. Normarily this would have raised an error, but as the type of 's' is 'int' in your example, and thus bitwise XOR is a well defined operation (though not the one you want). The second mistake is actually my fault, as the code on the mailing list that you linked to was actually written by me and contains a few typos. I realize now that we can simplify this code and make it more general: def in_sheet(site): return site.pos[1] == 0 # site in x-z plane, y position == 0 sheet = [idx for idx, site in enumerate(sysf.sites) if in_sheet(site)] # remember that 'wf(0)' returns a **sequence of wavefunctions** # so we first need to index which wavefunction we want (here the 1st one) wavefunction[0][sheet] Here we first construct a *list* of indices of all the sites in the x-z plane, and then index the wavefunction (a numpy array) with this list of indices. The second part, i.e. plotting the returned wavefunction, requires using some features of Kwant that, while documented, are rarely used directly by end-users. TL;DR the following code snippet will give you what you want: # remember that 'wf(0)' returns a **sequence of wavefunctions** # so we first need to index which wavefunction we want (here the 1st one) wavefunction_in_plane = wavefunction[0][sheet] density = np.abs(wavefunction_in_plane)**2 positions = [site.pos for site in sysf.sites if in_sheet(site)] # map positions to x-z plane positions_xz = np.array([(x, z) for x, y, z in positions]) data, min, max = kwant.plotter.mask_interpolate(positions_xy, density) image = pyplot.imshow(data, origin='lower') pyplot.colorbar(image) pyplot.show() We use the function 'mask_interpolate' from 'kwant.plotter' that will return data that can be directly shown as a colormap with 'pyplot.imshow'. We also need to map the site positions to the x-z plane, as 'imshow' will not know how to plot the output of 'mask_interpolate' if we provide positions in R^3 to this latter function. I have also attached a working example. Happy Kwanting, Joe
Dear John, The solution I am proposing is not very clean, but at least it works fine. To plot the wavefunction of a system for a fixed slice, you need to save the result of the wavefunction for the desired sites, and then define a new 2D system recovering the slice you chose, and by reordering the results of the wavefunction to fit the new order of the sites in the new system, the plot becomes straightforward. for your program it works like this: # plots for the mode 0 and 1 sys = make_system_cyl(a=1, s=50, L = 15) kwant.plot(sys) sysf = sys.finalized() params = SimpleNamespace(t=1) wf = kwant.wave_function(sysf, 1.5, args=[params]) #define a new system for the desired slice. lat2 = kwant.lattice.square() sys_b = kwant.Builder() #new 2D system used for the plot slice_x=5 mode=1 def Plot_mode(sys_b,slice_x,mode): wavefunction = wf(1)[mode] Sites= [site for (index,site) in enumerate(sysf.sites) if site.pos[0]==slice_x ] Indexes=[index for (index,site) in enumerate(sysf.sites) if site.pos[0]==slice_x ] waves=[wavefunction[index] for index in Indexes] Sites2=[lat2(site.pos[1],site.pos[2]) for site in Sites] #used to define the new system sys_b sys_b[Sites2]=0 sys_b[lat2.neighbors()]=-1 Sites_b=[site for site in sys_b.finalized().sites] # sites of the new system: with a different order! Index_b=[Sites2.index(site) for site in Sites_b] waves_b=[waves[index] for index in Index_b] # results is reordered according to the new order of the sites. return sys_b, waves_b sys_b, waves_b= Plot_mode(sys_b=sys_b,slice_x=5,mode=0) kwant.plotter.map(sys_b.finalized(),abs(array(waves_b))**2) sys_b, waves_b= Plot_mode(sys_b=sys_b,slice_x=5,mode=1) kwant.plotter.map(sys_b.finalized(),abs(array(waves_b))**2) I hope this helps. Regards, Adel On Tue, Apr 25, 2017 at 5:43 PM, John Eaves <johneaves02@gmail.com> wrote:
Good day,
I'm trying to reconstruct something discussed in http://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00432.html. In this answer it is described how to take the wavefunction from a 3D system and plot it in a 2D slice. This turned out to be something different than what the OP was asking for, but that is not so important. It is also mentioned that this has been asked before, but I couldn't find the right search terms. If my question has already been answered, I'd be very grateful for a referral to the relevant discussion!
Now, as it turns out, I can't seem to get the above to work. I take some small cylinder, I calculate the wavefunction, but then I get stuck at how to extract for example the y = 0 plane and plot this in 2D. Below I'll post how far I got, perhaps the fix would be rather simple.
def onsite(site, p): return 6*p.t
def hopping(site0, site1, p): return -p.t
def make_system_cyl(a=1, s=10, L = 15): lat = kwant.lattice.general([(0, 0, a), (a, 0, 0), (0, a, 0)])
sys = kwant.Builder()
def cyl_shape(pos): x, y, z = pos rsq = y ** 2 + z ** 2 return rsq < s^2 and 0 <= x < L
sys[lat.shape(cyl_shape, (1,1,1))] = onsite sys[lat.neighbors()] = hopping
def lead_shape(pos): x, y, z = pos rsq = y ** 2 + z ** 2 return rsq < s^2
sym_lead = kwant.TranslationalSymmetry((-a,0,0)) lead = kwant.Builder(sym_lead)
lead[lat.shape(lead_shape, (-a,1,1))] = onsite lead[lat.neighbors()] = hopping
sys.attach_lead(lead) sys.attach_lead(lead.reversed())
return sys
and I extract the wavefunction at some energy
sys = make_system_cyl(a=1, s=10, L = 15) kwant.plot(sys) sysf = sys.finalized() params = SimpleNamespace(t=1) wf = kwant.wave_function(sysf, 1.5, args=[params])
And then I try to extract the relevant parts
wavefunction = wf(0) projected_wavefunction = np.empty((wavefunction.size, 1), dtype=np.complex)
def in_sheet(index, site): return site.pos[1] == 0 # site in x-z plane
sheet_idx_sites = filter(in_sheet, enumerate(sys.sites)
for i, (idx, site) in enumerate(sheet_idx_sites): projected_wavefunction[i] = wavefunction[idx]
Now this last part I took from the post referenced at the start, and it gives a syntax error. I'm not really sure how to fix it as it looks fine to me. What's wrong there?
And more importantly, once this is fixed, how can one then make a kwant.plotter.map type of plot at this x-z plane?
-- Abbout Adel
participants (3)
-
Abbout Adel -
John Eaves -
Joseph Weston