How to extract a sheet from 3D system?
Hello all, I want to know how I can extract a sheet from any 3D system, e.g. I have coded 3D BHZ model, found wavefunction inside the scattering region, I want to extract data for a sheet (e.g. y=1 surface) out of wavefunction array. But I do not know how each site is labeled. I mean, I want to know how sites are ordered or how each element in wavefunction array is related to each site of system. You can find my code here: http://nbviewer.ipython.org/urls/gist.githubusercontent.com/AliAsgharpour/99... Best regards, Ali
Hi, The finalized system has a `sites` attribute which is a sequence of sites ordered in the same way as in the wavefunction. You should know that in the wavefunction the values on degrees of freedom which belong to the same site are stored sequentially. I see you have 4 degrees of freedom per site in your model, therefore you can do something like the following: import numpy as np import kwant sys = make_system() fsys = sys.finalized() wavefunction = .... ... # get sites in "sheet" def in_sheet(index, site): return site.pos[1] == 0 # site in x-z plane sheet_idx_sites = filter(in_sheet, enumerate(fsys.sites)) n = len(sheet_sites) # prepare storage projected_wavefunction = np.empty((n, 4), dtype=np.complex) # get wavefunction projection for i, (idx, site) in enumerate(sheet_idx_sites): projected_wavefunction[i] = wavefunction[idx*4: (idx+1)*4] Similar questions have been asked several times on the mailing list, although never explicitly with >1 orbital per site. Hope that helps, Joe P.S. I'm actually looking to implement better support for manipulating quantities which are defined over sites/hoppings, and am talking with Christoph about having this in kwant 2.0.
Thanks a lot Joseph. As you found my system has four degrees of freedom per site and I know (by reading previous emails) they are located sequentially in wavefunction array. I wanted to know how to filter them to access specific sites. Thanks again for your help, but I have another question about what you wrote. What is index in mentioned function "in_sheet"? Would you please explain how I can implement it in my code or a simple system? Thanks for your time and consideration. Best regards, Ali On Wed, Jul 29, 2015 at 8:06 PM, Joseph Weston <joseph.weston@cea.fr> wrote:
Hi,
The finalized system has a `sites` attribute which is a sequence of sites ordered in the same way as in the wavefunction. You should know that in the wavefunction the values on degrees of freedom which belong to the same site are stored sequentially.
I see you have 4 degrees of freedom per site in your model, therefore you can do something like the following:
import numpy as np import kwant
sys = make_system() fsys = sys.finalized()
wavefunction = ....
...
# get sites in "sheet"
def in_sheet(index, site): return site.pos[1] == 0 # site in x-z plane
sheet_idx_sites = filter(in_sheet, enumerate(fsys.sites)) n = len(sheet_sites)
# prepare storage projected_wavefunction = np.empty((n, 4), dtype=np.complex)
# get wavefunction projection for i, (idx, site) in enumerate(sheet_idx_sites): projected_wavefunction[i] = wavefunction[idx*4: (idx+1)*4]
Similar questions have been asked several times on the mailing list, although never explicitly with >1 orbital per site.
Hope that helps,
Joe
P.S. I'm actually looking to implement better support for manipulating quantities which are defined over sites/hoppings, and am talking with Christoph about having this in kwant 2.0.
Hi, I think I misunderstood what you are trying to do. You talked about "filtering" only the sites which were in a plane, so I thought you wanted to produce a projection of the wavefunction onto the sites in the plane. If you just want to be able to access wavefunction elements "by their site" you could do something like the following def site_indexer(wavefunction, sites, n_orbs=4): site_idx_map = {site: i for i, site in enumerate(sites)} def indexer(site): orbital = site_idx_map[site] * n_orbs return wavefunction[orbital: orbital + n_orbs] return indexer lat = kwant.lattice.square() sys = make_system() fsys = sys.finalized() wf = kwant.wavefunction(...) better_wf = site_indexer(wf, fsys.sites) # we get an array of length 4, for the 4 orbitals on the site print better_wf(lat(0, 2)) You could then use `better_wf` to access the wavefunction on the orbitals belonging to different sites. Joe
Thanks again for your prompt answer. I want to access wavefunction elements by their site which are located in a plane for 3D or a line in 2D, I mean no projected wavefunction, just those elements which are related to sites on a plane in 3D system or a line in 2D one. I implemented what you said for 2D system, but I don't get any array for "better_wf(lat(i.j))". It should give me an array (4 elements due to 4 degrees of freedom) for a specific site. If it works, I can change it to the ones I want, but it does not. Here you can find my code: http://nbviewer.ipython.org/urls/gist.githubusercontent.com/AliAsgharpour/2d... Your time and consideration are much appreciated. Best regards, Ali On Thu, Jul 30, 2015 at 5:04 PM, Joseph Weston <joseph.weston@cea.fr> wrote:
Hi,
I think I misunderstood what you are trying to do. You talked about "filtering" only the sites which were in a plane, so I thought you wanted to produce a projection of the wavefunction onto the sites in the plane. If you just want to be able to access wavefunction elements "by their site" you could do something like the following
def site_indexer(wavefunction, sites, n_orbs=4): site_idx_map = {site: i for i, site in enumerate(sites)}
def indexer(site): orbital = site_idx_map[site] * n_orbs return wavefunction[orbital: orbital + n_orbs]
return indexer
lat = kwant.lattice.square() sys = make_system() fsys = sys.finalized() wf = kwant.wavefunction(...)
better_wf = site_indexer(wf, fsys.sites)
# we get an array of length 4, for the 4 orbitals on the site print better_wf(lat(0, 2))
You could then use `better_wf` to access the wavefunction on the orbitals belonging to different sites.
Joe
Yeah, ok so `kwant.wavefunction(...)` returns a callable which you call with a lead and you get back an array which we will call `wfs`. `wfs` is a *2D* array, the first axis is the *mode number*, the second axis is the *orbital number*. See the docs: http://kwant-project.org/doc/1.0/reference/generated/kwant.solvers.default.w... So you need to do: wf = kwant.wave_function(fsys, energy=1.7, args=[par,]) wf0, wf1 = wavefunctions(0) # iterate over first (mode) axis better_wf0 = site_indexer(wf0) better_wf1 = site_indexer(wf1) better_wf0(lat(1,0)) # returns a 4-element array I checked it with your notebook and it works Regards, Joe
Thank you so much. I really appreciate your help. Best regards, Ali On Thu, Jul 30, 2015 at 7:05 PM, Joseph Weston <joseph.weston@cea.fr> wrote:
Yeah, ok so `kwant.wavefunction(...)` returns a callable which you call with a lead and you get back an array which we will call `wfs`. `wfs` is a *2D* array, the first axis is the *mode number*, the second axis is the *orbital number*. See the docs:
http://kwant-project.org/doc/1.0/reference/generated/kwant.solvers.default.w...
So you need to do:
wf = kwant.wave_function(fsys, energy=1.7, args=[par,]) wf0, wf1 = wavefunctions(0) # iterate over first (mode) axis
better_wf0 = site_indexer(wf0) better_wf1 = site_indexer(wf1)
better_wf0(lat(1,0)) # returns a 4-element array
I checked it with your notebook and it works
Regards,
Joe
participants (2)
-
Ali Asgharpour -
Joseph Weston