Visualize current through a 2D cut for a 3D system
Hello everyone, Is there currently way to visualize the current through a 2D cut for a 3D system that does not involve messing with the kwant code? I have a 3D system and would like to see the differences in the current when I switch some specific hoppings on and off. Regards,  Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
Dear Eleni, I think this should be doable with not too much of low level work. * Firstly you probably will want to select only the hoppings that are close to the 2D cut. You can do this by using "where" parameter to Current operator (see https://kwantproject.org/doc/1/reference/generated/kwant.operator.Current#k...) * After computing the current across those hoppings, you would need to call kwant.plotter.interpolate_current. This function only does the interpolation, but no plotting. Having that interpolation in 3D you need to slice it and probably take a projection onto 2D plane. * Finally, the resulting 2D array you can feed to kwant.plotter.streamplot kwant.plotter.current is essentially a straightforward combination of kwant.plotter.interpolate_current and kwant.plotter.streamplot Also: your problem sounds like a useful thing to do. Please let me know if you succeed, and what you ended up doing (or if you face any further problems as well). Best, Anton On Fri, Mar 23, 2018 at 5:22 PM, <elchatz@auth.gr> wrote:
Hello everyone,
Is there currently way to visualize the current through a 2D cut for a 3D system that does not involve messing with the kwant code?
I have a 3D system and would like to see the differences in the current when I switch some specific hoppings on and off.
Regards,
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
Hello Anton, Thank you very much for your reply. I am having a hard time with the "where" argument in kwant.operator.Current. My system is not finalized. How can I find the sites associated with the hopping? By their tag? Their position? Do you have a piece of code for that? Note that I am using TBModels and Wannier90 and the lattice has one sublattice for each WF (or site, or orbital). norb is not defined, but I see that for x in lattice.sublattices: x.norbs=1 works and looks correct. Regards, Eleni Quoting Anton Akhmerov <anton.akhmerov+kd@gmail.com>:
Dear Eleni,
I think this should be doable with not too much of low level work. * Firstly you probably will want to select only the hoppings that are close to the 2D cut. You can do this by using "where" parameter to Current operator (see https://kwantproject.org/doc/1/reference/generated/kwant.operator.Current#k...) * After computing the current across those hoppings, you would need to call kwant.plotter.interpolate_current. This function only does the interpolation, but no plotting. Having that interpolation in 3D you need to slice it and probably take a projection onto 2D plane. * Finally, the resulting 2D array you can feed to kwant.plotter.streamplot
kwant.plotter.current is essentially a straightforward combination of kwant.plotter.interpolate_current and kwant.plotter.streamplot
Also: your problem sounds like a useful thing to do. Please let me know if you succeed, and what you ended up doing (or if you face any further problems as well).
Best, Anton
On Fri, Mar 23, 2018 at 5:22 PM, <elchatz@auth.gr> wrote:
Hello everyone,
Is there currently way to visualize the current through a 2D cut for a 3D system that does not involve messing with the kwant code?
I have a 3D system and would like to see the differences in the current when I switch some specific hoppings on and off.
Regards,
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
Anton Akhmerov wrote:
On Fri, Mar 23, 2018 at 5:22 PM, <elchatz@auth.gr> wrote:
Is there currently way to visualize the current through a 2D cut for a 3D system that does not involve messing with the kwant code?
* After computing the current across those hoppings, you would need to call kwant.plotter.interpolate_current. This function only does the interpolation, but no plotting. Having that interpolation in 3D you need to slice it and probably take a projection onto 2D plane.
This gives me the following idea: Currently, interpolate_current interpolates the current in a ndimensional Kwant system onto a ndimensional grid. But as we can from Eleni's request it can be useful to generalize this. Without complicating the interface for simple use cases we could do the following: we separate the interpolation algorithm from interpolate_current. I.e. there would be a function (e.g. interpolate_vector_field) that takes a ndimensional current field (i.e. a sequence of pairs of points together with the current strengths) and interpolates this onto a given mdimensional grid, with m <= n. Without loss of generality, the grid can be oriented along the first m axes since an arbitrary transformation may be applied to the hoppings before interpolate_vector_field is called. This would allow Eleni to plot currents across any plane. Incidentally, the function interpolate_vector_field matches quite exactly the part of interpolate_current that should be implemented in C to speed it up.
Hello everyone, I'm not sure what I am doing wrong, but here is the output when using interpolate_current(): 
kwant_sys = wraparound.wraparound(kwant_model).finalized() ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) evecs = sla.eigsh(ham_mat, k=48, which='SM')[1] J = kwant.operator.Current(kwant_sys) current = J(evecs[:, 9]) IJ = kwant.plotter.interpolate_current(kwant_sys, J) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/apps/applications/python/anaconda3/5.0.1/lib/python3.6/sitepackages/kwant/plotter.py", line 1852, in interpolate_current if len(current) != syst.graph.num_edges: TypeError: object of type 'kwant.operator.Current' has no len() current array([ 1.04495760e10, 1.04407946e10, 4.08360326e10, ..., 5.89810275e10, 2.77384947e10, 1.41194501e10])
I was also looking at this thread: https://mailmanmail5.webfaction.com/pipermail/kwantdiscuss/2017October/00... and I think I will try also using ipyvolume or mayavi. But I would like to know if there is a problem with this being a closed system. Regards, Eleni Quoting Anton Akhmerov <anton.akhmerov+kd@gmail.com>:
Dear Eleni,
I think this should be doable with not too much of low level work. * Firstly you probably will want to select only the hoppings that are close to the 2D cut. You can do this by using "where" parameter to Current operator (see https://kwantproject.org/doc/1/reference/generated/kwant.operator.Current#k...) * After computing the current across those hoppings, you would need to call kwant.plotter.interpolate_current. This function only does the interpolation, but no plotting. Having that interpolation in 3D you need to slice it and probably take a projection onto 2D plane. * Finally, the resulting 2D array you can feed to kwant.plotter.streamplot
kwant.plotter.current is essentially a straightforward combination of kwant.plotter.interpolate_current and kwant.plotter.streamplot
Also: your problem sounds like a useful thing to do. Please let me know if you succeed, and what you ended up doing (or if you face any further problems as well).
Best, Anton
On Fri, Mar 23, 2018 at 5:22 PM, <elchatz@auth.gr> wrote:
Hello everyone,
Is there currently way to visualize the current through a 2D cut for a 3D system that does not involve messing with the kwant code?
I have a 3D system and would like to see the differences in the current when I switch some specific hoppings on and off.
Regards,
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
Hi, On 04/10/2018 01:50 PM, elchatz@auth.gr wrote:
Hello everyone,
I'm not sure what I am doing wrong, but here is the output when using interpolate_current():

kwant_sys = wraparound.wraparound(kwant_model).finalized() ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) evecs = sla.eigsh(ham_mat, k=48, which='SM')[1] J = kwant.operator.Current(kwant_sys) current = J(evecs[:, 9]) IJ = kwant.plotter.interpolate_current(kwant_sys, J)
You're passing "J" to "interpolate current", when you should pass "current"
and I think I will try also using ipyvolume or mayavi. But I would like to know if there is a problem with this being a closed system.
No, no problem. Of course if you don't break time reversal symmetry (e.g. with magnetic field) then you're eigenstates will have 0 current everywhere. Happy Kwanting, Joe
Hello Joseph, Sorry I copied the wrong pieces of code:  ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) ev = sla.eigsh(ham_mat, k=48, which='SM') evecs = ev[1] lat1 = lattice.sublattices[1] lat3 = lattice.sublattices[3] J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0, 0, 0))]) current = J(evecs[:, 1]) IJ = kwant.plotter.interpolate_current(kwant_sys, current)  ValueError Traceback (most recent call last) <ipythoninput162090a8a157546> in <module>() 5 J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0, 0, 0))]) 6 current = J(evecs[:, 1]) > 7 IJ = kwant.plotter.interpolate_current(kwant_sys, current) 8 print(current) 9 ~\Anaconda3\lib\sitepackages\kwant\plotter.py in interpolate_current(syst, current, relwidth, abswidth, n) 1851 1852 if len(current) != syst.graph.num_edges: > 1853 raise ValueError("Current and hoppings arrays do not have the same" 1854 " length.") 1855 ValueError: Current and hoppings arrays do not have the same length. print(current) [7.65418274e10]   ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) ev = sla.eigsh(ham_mat, k=48, which='SM') evecs = ev[1] J = kwant.operator.Current(kwant_sys) current = J(evecs[:, 1]) IJ = kwant.plotter.interpolate_current(kwant_sys, current)  ValueError Traceback (most recent call last) <ipythoninput164ae3dc103f55d> in <module>() 6 J = kwant.operator.Current(kwant_sys) 7 current = J(evecs[:, 1]) > 8 IJ = kwant.plotter.interpolate_current(kwant_sys, current) 9 ~\Anaconda3\lib\sitepackages\kwant\plotter.py in interpolate_current(syst, current, relwidth, abswidth, n) 1900 # this check. This check is done here to keep changes local 1901 if dim != 2: > 1902 raise ValueError("'interpolate_current' only works for 2D systems.") 1903 factor = (3 / np.pi) / (width / 2) 1904 scale = 2 / width ValueError: 'interpolate_current' only works for 2D systems.  I was able to do a points3d() plot with the 3D data in Mayavi using the following piece of code:  ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) ev = sla.eigsh(ham_mat, k=48, which='SM') evecs = ev[1] prob_dens = np.abs(evecs[:, 47])**2 x = np.array([]) y = np.array([]) z = np.array([]) for i in range(len(kwant_sys.sites)): x = np.append(x,[kwant_sys.pos(i)[0]], axis=0) y = np.append(y,[kwant_sys.pos(i)[1]], axis=0) z = np.append(z,[kwant_sys.pos(i)[2]], axis=0) pts = mlab.points3d(x, y, z, prob_dens) mayavi.mlab.scalarbar(object=pts) mayavi.mlab.xlabel("X", object = pts) mlab.savefig("pts_47.png") mlab.show()  However, a contour3d would be more appropriate, but it looks like that requires regularily spaced points. On the other hand, an unstructured grid looks a bit more cumbersome to create. At this point, the easiest thing looks like copying from kwant code directly into my notebook for interpolation etc. Eleni Quoting Joseph Weston <joseph.weston08@gmail.com>:
Hi,
On 04/10/2018 01:50 PM, elchatz@auth.gr wrote:
Hello everyone,
I'm not sure what I am doing wrong, but here is the output when using interpolate_current():

kwant_sys = wraparound.wraparound(kwant_model).finalized() ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) evecs = sla.eigsh(ham_mat, k=48, which='SM')[1] J = kwant.operator.Current(kwant_sys) current = J(evecs[:, 9]) IJ = kwant.plotter.interpolate_current(kwant_sys, J)
You're passing "J" to "interpolate current", when you should pass "current"
and I think I will try also using ipyvolume or mayavi. But I would like to know if there is a problem with this being a closed system.
No, no problem. Of course if you don't break time reversal symmetry (e.g. with magnetic field) then you're eigenstates will have 0 current everywhere.
Happy Kwanting,
Joe
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
Hi,
Hello Joseph,
Sorry I copied the wrong pieces of code:
 ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) ev = sla.eigsh(ham_mat, k=48, which='SM') evecs = ev[1] lat1 = lattice.sublattices[1] lat3 = lattice.sublattices[3] J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0, 0, 0))]) current = J(evecs[:, 1]) IJ = kwant.plotter.interpolate_current(kwant_sys, current) 
ValueError Traceback (most recent call last) <ipythoninput162090a8a157546> in <module>() 5 J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0, 0, 0))]) 6 current = J(evecs[:, 1]) > 7 IJ = kwant.plotter.interpolate_current(kwant_sys, current) 8 print(current) 9
~\Anaconda3\lib\sitepackages\kwant\plotter.py in interpolate_current(syst, current, relwidth, abswidth, n) 1851 1852 if len(current) != syst.graph.num_edges: > 1853 raise ValueError("Current and hoppings arrays do not have the same" 1854 " length.") 1855
ValueError: Current and hoppings arrays do not have the same length.
print(current) [7.65418274e10] 
Ah right, you are only calculating the current over a single hopping, but 'interpolate_current' requires the current over *all* hoppings. What extra information do you hope to get by plotting in realspace the current across a single hopping? This reason this restriction is in place is that 'interpolate_current' accepts just an array of numbers for the current. This means that there is no way for 'interpolate_current' to know which number corresponds to the current across which hopping without further information. This 'further information' is the fact that the array of currents is in the same order as 'syst.graph'; if we don't specify all the hoppings there's know way to know which ones we've specified.
 ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) ev = sla.eigsh(ham_mat, k=48, which='SM') evecs = ev[1] J = kwant.operator.Current(kwant_sys) current = J(evecs[:, 1]) IJ = kwant.plotter.interpolate_current(kwant_sys, current) 
ValueError Traceback (most recent call last) <ipythoninput164ae3dc103f55d> in <module>() 6 J = kwant.operator.Current(kwant_sys) 7 current = J(evecs[:, 1]) > 8 IJ = kwant.plotter.interpolate_current(kwant_sys, current) 9
~\Anaconda3\lib\sitepackages\kwant\plotter.py in interpolate_current(syst, current, relwidth, abswidth, n) 1900 # this check. This check is done here to keep changes local 1901 if dim != 2: > 1902 raise ValueError("'interpolate_current' only works for 2D systems.") 1903 factor = (3 / np.pi) / (width / 2) 1904 scale = 2 / width
ValueError: 'interpolate_current' only works for 2D systems.

The restriction of 'interpolate_current' to 2D was lifted a while ago (https://gitlab.kwantproject.org/kwant/kwant/merge_requests/166) but we have not made an official release since. If you're feeling adventurous you can install the Kwant directly from the 'master' branch of the git repository (post back here if you don't know how). Happy Kwanting, Joe
Hello all, Is it also possible to help me with this function to match first nearest neighbor sites for the where argument in kwant.operator.Current: def first_neighbors(site_to, site_from): for i in site_from.family.neighbors(1): if(??): return True else return False [...] J = kwant.operator.Current(kwant_sys, where=first_neighbors, sum=True] Thank you Quoting Joseph Weston <joseph.weston08@gmail.com>:
Hi,
Hello Joseph,
Sorry I copied the wrong pieces of code:
 ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) ev = sla.eigsh(ham_mat, k=48, which='SM') evecs = ev[1] lat1 = lattice.sublattices[1] lat3 = lattice.sublattices[3] J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0, 0, 0))]) current = J(evecs[:, 1]) IJ = kwant.plotter.interpolate_current(kwant_sys, current) 
ValueError Traceback (most recent call last) <ipythoninput162090a8a157546> in <module>() 5 J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0, 0, 0))]) 6 current = J(evecs[:, 1]) > 7 IJ = kwant.plotter.interpolate_current(kwant_sys, current) 8 print(current) 9
~\Anaconda3\lib\sitepackages\kwant\plotter.py in interpolate_current(syst, current, relwidth, abswidth, n) 1851 1852 if len(current) != syst.graph.num_edges: > 1853 raise ValueError("Current and hoppings arrays do not have the same" 1854 " length.") 1855
ValueError: Current and hoppings arrays do not have the same length.
print(current) [7.65418274e10] 
Ah right, you are only calculating the current over a single hopping, but 'interpolate_current' requires the current over *all* hoppings. What extra information do you hope to get by plotting in realspace the current across a single hopping?
This reason this restriction is in place is that 'interpolate_current' accepts just an array of numbers for the current. This means that there is no way for 'interpolate_current' to know which number corresponds to the current across which hopping without further information. This 'further information' is the fact that the array of currents is in the same order as 'syst.graph'; if we don't specify all the hoppings there's know way to know which ones we've specified.
 ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) ev = sla.eigsh(ham_mat, k=48, which='SM') evecs = ev[1] J = kwant.operator.Current(kwant_sys) current = J(evecs[:, 1]) IJ = kwant.plotter.interpolate_current(kwant_sys, current) 
ValueError Traceback (most recent call last) <ipythoninput164ae3dc103f55d> in <module>() 6 J = kwant.operator.Current(kwant_sys) 7 current = J(evecs[:, 1]) > 8 IJ = kwant.plotter.interpolate_current(kwant_sys, current) 9
~\Anaconda3\lib\sitepackages\kwant\plotter.py in interpolate_current(syst, current, relwidth, abswidth, n) 1900 # this check. This check is done here to keep changes local 1901 if dim != 2: > 1902 raise ValueError("'interpolate_current' only works for 2D systems.") 1903 factor = (3 / np.pi) / (width / 2) 1904 scale = 2 / width
ValueError: 'interpolate_current' only works for 2D systems.

The restriction of 'interpolate_current' to 2D was lifted a while ago (https://gitlab.kwantproject.org/kwant/kwant/merge_requests/166) but we have not made an official release since. If you're feeling adventurous you can install the Kwant directly from the 'master' branch of the git repository (post back here if you don't know how).
Happy Kwanting,
Joe
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
Hi,
Is it also possible to help me with this function to match first nearest neighbor sites for the where argument in kwant.operator.Current:
def first_neighbors(site_to, site_from): for i in site_from.family.neighbors(1): if(??): return True else return False
[...]
J = kwant.operator.Current(kwant_sys, where=first_neighbors, sum=True]
Probably the following: def first_neighbors(to_site, from_site): return to_site in from_site.neighbors(1) 'neighbors' returns (per the documentation [1]) a list of sites and we can just use the 'in' operator to check for the requested site. Happy Kwanting, Joe [1]: https://kwantproject.org/doc/1.0/reference/generated/kwant.lattice.Polyatom...
Hello Joseph, One problem that I see is that I can only get neighbors() from Builders and lattices (using site.family). Using something like this: def first_neighbors(to_site, from_site): return to_site in from_site.family.neighbors(1) is a good idea, but am I really comparing two sites?  for i in x.family.neighbors(1): print(i) HoppingKind((1, 0, 0), <Monatomic lattice 13>) HoppingKind((0, 1, 0), <Monatomic lattice 13>)  Also, since this is very expensive computationally, is there any other way to find the hopping distance between two sites? Eleni Quoting Joseph Weston <joseph.weston08@gmail.com>:
Hi,
Is it also possible to help me with this function to match first nearest neighbor sites for the where argument in kwant.operator.Current:
def first_neighbors(site_to, site_from): for i in site_from.family.neighbors(1): if(??): return True else return False
[...]
J = kwant.operator.Current(kwant_sys, where=first_neighbors, sum=True]
Probably the following:
def first_neighbors(to_site, from_site): return to_site in from_site.neighbors(1)
'neighbors' returns (per the documentation [1]) a list of sites and we can just use the 'in' operator to check for the requested site.
Happy Kwanting,
Joe
[1]: https://kwantproject.org/doc/1.0/reference/generated/kwant.lattice.Polyatom...
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
Hi again,
Hello Joseph,
One problem that I see is that I can only get neighbors() from Builders and lattices (using site.family).
Using something like this:
def first_neighbors(to_site, from_site): return to_site in from_site.family.neighbors(1)
is a good idea, but am I really comparing two sites?
 for i in x.family.neighbors(1): print(i)
HoppingKind((1, 0, 0), <Monatomic lattice 13>) HoppingKind((0, 1, 0), <Monatomic lattice 13>) 
Ah, sorry I was confused.
Also, since this is very expensive computationally, is there any other way to find the hopping distance between two sites?
Sure! If the sites are from the same lattice then you can just take the difference of the "tags" of the sites: delta = to_site.tag  from_site.tag 'delta' will be an integer array: the position of the site in the basis of lattice vectors. If the sites are from different lattices then you probably can't do better than comparing the realspace positions: realspace_delta = to_site.pos  from_site.pos Happy Kwanting, Joe
Thank you. That's great. Quoting Joseph Weston <joseph.weston08@gmail.com>:
Hi again,
Hello Joseph,
One problem that I see is that I can only get neighbors() from Builders and lattices (using site.family).
Using something like this:
def first_neighbors(to_site, from_site): return to_site in from_site.family.neighbors(1)
is a good idea, but am I really comparing two sites?
 for i in x.family.neighbors(1): print(i)
HoppingKind((1, 0, 0), <Monatomic lattice 13>) HoppingKind((0, 1, 0), <Monatomic lattice 13>) 
Ah, sorry I was confused.
Also, since this is very expensive computationally, is there any other way to find the hopping distance between two sites?
Sure! If the sites are from the same lattice then you can just take the difference of the "tags" of the sites:
delta = to_site.tag  from_site.tag
'delta' will be an integer array: the position of the site in the basis of lattice vectors. If the sites are from different lattices then you probably can't do better than comparing the realspace positions:
realspace_delta = to_site.pos  from_site.pos
Happy Kwanting,
Joe
 Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elchatz@auth.gr  tel:+30 2310 998109
participants (4)

Anton Akhmerov

Christoph Groth

elchatz＠auth.gr

Joseph Weston