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 <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