# Re: [Kwant] Visualize current through a 2D cut for a 3D system

```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.westo...@gmail.com>:

```
```Hi,

On 04/10/2018 01:50 PM, elch...@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
elch...@auth.gr - tel:+30 2310 998109

```