Draw Bond-Currents with sisl

[Mr-Bey]

Hi Nick,

I got a tutorial on calculating the bond-current from your answer here. And I tried to draw the bonding current of graphene. But the atom is not shown in the graph, and the image does not seem to be correct.

This is the code:

import os
import numpy as np
import sisl

# Now read in the AV file
av = sisl.get_sile('graphene_scat.TBT.AV.nc')
geom = av.geom

J_orb = av.orbital_current('Left', E=2, isc=[None] * 3)
J_a = av.atom_current_from_orbital(J_orb, activity=False)

#J_v = av.vector_current_from_orbital(J_orb)

# Now create plot
from matplotlib.colors import Colormap
from matplotlib.mlab import griddata
import matplotlib.pyplot as plt

def mmlinspace(a, axis, N):
    m, M = np.min(a[:, axis]), np.max(a[:, axis])
    return np.linspace(m, M, N)

N_x, N_y = 2000, 2000
x = mmlinspace(geom.xyz, 0, N_x)
y = mmlinspace(geom.xyz, 1, N_y)

xi, yi = np.meshgrid(x, y)
zi = griddata(geom.xyz[:, 0], geom.xyz[:, 1], J_a, xi, yi, interp='linear')

cm = Colormap('inferno')
plt.contourf(xi, yi, zi, cmap=plt.cm.plasma)
plt.savefig('test2.png')

plt.show()

My question is how to set "J_orb" correct in the code, and If I expect to get both the atom and bond-current, What changes should I make for this code?

Best Regards~

[zerothi]

What you are plotting in that routine is the atomic currents (not bond-currents).

The bond-currents may be retrieved by:

av = sisl.get_sile('graphene_scat.TBT.AV.nc')
geom = av.geom
J_bond = av.bond_current('Left', E=2.0)

Now J_bond is the sparse matrix with bond-currents that you need to plot.

So basically you need to plot arrows between atoms now. You should decide on how to plot this.

In https://doi.org/10.1088/1361-648X/aad6f1 there is a section describing how to plot the bond-currents.

[Mr-Bey]

When using this code:

geom = av.geom
J_orb = av.bond_current('Left', E=1.98)
J_a = av.bond_current_from_orbital(J_orb)

I get the following error message:

Traceback (most recent call last):
File "read_tbt.py", line 11, in <module>
J_a = av.bond_current_from_orbital(J_orb)
File "/home/lee/anaconda3/lib/python3.6/site-packages/sisl/io/tbtrans/tbt.py", line 1150, in bond_current_from_orbital 
indptr = np.insert(_a.cumsumi(iptr[fo[1:]] - iptr[fo[:-1]]), 0, 0)
IndexError: index 168 is out of bounds for axis 0 with size 161

[zerothi]

J_orb from the above code is actually the bond-current. So no need to pass it through to the next call.

Please carefully read the documentation of: atom_current, bond_current and orbital_current in the documentation of the tbtrans output file: http://zerothi.github.io/sisl/api-generated/sisl.io.tbtrans.tbtncSileTBtrans.html#sisl.io.tbtrans.tbtncSileTBtrans

[Mr-Bey]

I'm not good at the tool matplotlib, so I exported the sparse array and geometric data， then used Matlab to draw the picture. But the picture drawn is completely different from the one in the literature (the structure of graphene is later superimposed).

Would you mind give me some help in the code about drawing this sparse array?

[zerothi]

Generally I think you should figure out how to do this.

I can provide this small code which plots the bond-currents. However, note that it is very ugly so you should refine and adapt for your needs.

import numpy as np
import sisl as s
import matplotlib.pyplot as plt

# Get TBT file
tbt = s.get_sile('siesta.TBT.nc')
geom = tbt.geometry

# Retrieve bond-currents
J = tbt.bond_current(0, 0.1)

# Get X and Y coordinates
x, y = geom.xyz[:, :2].T

# Create plot
fig = plt.figure()
ax = fig.add_subplot(111)

# Add geometry
ax.scatter(x, y)

def r(N, x):
    return np.repeat(x, N)

# Add quiver plots
scale = 0.1
for ia in geom:
    Jia = J[ia, :]
    if Jia.nnz == 0:
        # There are no bond-currents here
        continue
    # Calculate vectors
    Rx = geom.axyz(Jia.indices)[:, 0] - x[ia]
    Ry = geom.axyz(Jia.indices)[:, 1] - y[ia]
    d = (Rx ** 2 + Ry ** 2) ** 0.5
    Rx *= Jia.data[:] / d * scale
    Ry *= Jia.data[:] / d * scale

    N = len(Rx)
    ax.quiver(r(N, x[ia]), r(N, y[ia]), Rx, Ry, Jia.data[:] / Jia.data.max())

plt.show()

[zerothi]

I will consider this issue resolved :)

[Mr-Bey]

Yes，it resolved. And thanks a lot for your help！Cheers！XD