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 :

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://kwant-project.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, 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

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/site-packages/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.04495760e-10, 1.04407946e-10, 4.08360326e-10, ..., -5.89810275e-10, -2.77384947e-10, 1.41194501e-10])

I was also looking at this thread: https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-October/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 :

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://kwant-project.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, 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

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) <ipython-input-162-090a8a157546> 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\site-packages\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.65418274e-10] -------------------------------------------------- --------------------------------------------------- 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) <ipython-input-164-ae3dc103f55d> 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\site-packages\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 :

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

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 :

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) <ipython-input-162-090a8a157546> 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\site-packages\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.65418274e-10] --------------------------------------------------

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) <ipython-input-164-ae3dc103f55d> 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\site-packages\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.kwant-project.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

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), ) HoppingKind((0, 1, 0), ) ------------------------------------------- Also, since this is very expensive computationally, is there any other way to find the hopping distance between two sites? Eleni Quoting Joseph Weston :

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://kwant-project.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

Thank you. That's great. Quoting Joseph Weston :

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), ) HoppingKind((0, 1, 0), ) -------------------------------------------

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