In [ ]:
from __future__ import print_function
import sisl
import numpy as np
from matplotlib import animation
import matplotlib.pyplot as plt
from IPython.display import HTML
%matplotlib inline
In this example you will learn how to make use of the periodicity of the electrodes.
As seen in TB 4 the transmission calculation takes a considerable amount of time. In this example we will redo the same calculation, but speed it up (no approximations made).
A large computational effort is made on calculating the self-energies which basically is inverting, multiplying and adding matrices, roughly 10-20 times per $k$-point, per energy point, per electrode.
For some systems this is far more demanding than calculating the Green function for the system.
In systems where there is periodicity along the transverse semi-infinite direction (not along the transport direction) one can utilize Bloch's theorem to reduce the computational cost of calculating the self-energy. In ANY calculation if you have periodicity, please USE it.
In this example you should scour the tbtrans manual on how to enable Bloch's theorem, and once enabled it should be roughly 3 - 4 times as fast, something that is non-negligeble for large systems.
In [ ]:
graphene = sisl.geom.graphene(orthogonal=True)
Note the below lines are differing from the same lines in TB 4, i.e. we save the electrode electronic structure without extending it 25 times.
In [ ]:
H_elec = sisl.Hamiltonian(graphene)
H_elec.construct(([0.1, 1.43], [0., -2.7]))
H_elec.write('ELEC.nc')
See TB 2 for details on why we choose repeat/tile on the Hamiltonian object and not on the geometry, prior to construction.
In [ ]:
H = H_elec.repeat(25, axis=0).tile(15, axis=1)
H = H.remove(H.geom.close(H.geom.center(what='cell'), R=10.))
dangling = [ia for ia in H.geom.close(H.geom.center(what='cell'), R=14.)
if len(H.edges(ia)) < 2]
H = H.remove(dangling)
# + 1 to get fortran indices
edge = [ia + 1 for ia in H.geom.close(H.geom.center(what='cell'), R=14.)
if len(H.edges(ia)) < 3]
# Pretty-print the list of atoms
print(sisl.utils.list2str(edge))
H.geom.write('device.xyz')
H.write('DEVICE.nc')
Instead of analysing the same thing as in TB 4 you should perform the following actions to explore the available data-analysis capabilities of TBtrans.
HINT please copy as much as you like from example 04 to simplify the following tasks.
tbt.In the following we will concentrate on only looking at $\Gamma$-point related quantities. I.e. all quantities should only be plotted for this $k$-point.
To extract information for one or more subset of points you should look into the function
help(tbt.kindex)
which may be used to find a resulting $k$-point index in the result file.
Plot the transmission ($\Gamma$-point only). To extract a subset $k$-point you should read the documentation for the functions (hint: kavg is the keyword you are looking for).
H_elec.tile(25, axis=0).tile(15, axis=1) instead of repeat/tile. Which of repeat or tile are faster?
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
tbt = sisl.get_sile('siesta.TBT.nc')
# Easier manipulation of the geometry
geom = tbt.geometry
a_dev = tbt.a_dev # the indices where we have DOS
# Extract the DOS, per orbital (hence sum=False)
DOS = tbt.ADOS(0, sum=False)
# Normalize DOS for plotting (maximum size == 400)
# This array has *all* energy points and orbitals
DOS /= DOS.max() / 400
a_xyz = geom.xyz[a_dev, :2]
In [ ]:
%%capture
fig = plt.figure(figsize=(12,4));
ax = plt.axes();
scatter = ax.scatter(a_xyz[:, 0], a_xyz[:, 1], 1);
ax.set_xlabel(r'$x$ [Ang]'); ax.set_ylabel(r'$y$ [Ang]');
ax.axis('equal');
In [ ]:
# If this animation does not work, then don't spend time on it!
def animate(i):
ax.set_title('Energy {:.3f} eV'.format(tbt.E[i]));
scatter.set_sizes(DOS[i]);
return scatter,
anim = animation.FuncAnimation(fig, animate, frames=len(tbt.E), interval=100, repeat=False)
HTML(anim.to_html5_video())
.edges for the Hamiltonian.*.TBT.nc file