In [ ]:
from __future__ import print_function
import sisl
import numpy as np
import matplotlib.pyplot as plt
from functools import partial
%matplotlib inline
TBtrans is capable of calculating transport in $N\ge 1$ electrode systems. In this example we will explore a 4-terminal graphene GNR cross-bar (one zGNR, the other aGNR) system.
In [ ]:
graphene = sisl.geom.graphene(orthogonal=True)
In [ ]:
R = [0.1, 1.43]
hop = [0., -2.7]
Create the two electrodes in $x$ and $y$ directions. We will force the systems to be nano-ribbons, i.e. only periodic along the ribbon. In sisl there are two ways of accomplishing this.
The below code uses the first method.
Please see if you can change the creation of elec_x by adding vacuum.
HINT: Look at the documentation for the sisl.Geometry and search for vacuum. To know the orbital distance look up maxR in the geometry class as well.
In [ ]:
elec_y = graphene.tile(3, axis=0)
elec_y.set_nsc([1, 3, 1])
elec_y.write('elec_y.xyz')
elec_x = graphene.tile(5, axis=1)
elec_x.set_nsc([3, 1, 1])
elec_x.write('elec_x.xyz')
Subsequently we create the electronic structure.
In [ ]:
H_y = sisl.Hamiltonian(elec_y)
H_y.construct((R, hop))
H_y.write('ELEC_Y.nc')
H_x = sisl.Hamiltonian(elec_x)
H_x.construct((R, hop))
H_x.write('ELEC_X.nc')
Now we have created the electronic structure for the electrodes. All that is needed is the electronic structure of the device region, i.e. the crossing nano-ribbons.
In [ ]:
dev_y = elec_y.tile(30, axis=1)
dev_y = dev_y.translate( -dev_y.center(what='xyz') )
dev_x = elec_x.tile(18, axis=0)
dev_x = dev_x.translate( -dev_x.center(what='xyz') )
Remove any atoms that are duplicated, i.e. when we overlay these two geometries some atoms are the same.
In [ ]:
device = dev_y.add(dev_x)
device.set_nsc([1,1,1])
duplicates = []
for ia in dev_y:
idx = device.close(ia, 0.1)
if len(idx) > 1:
duplicates.append(idx[1])
device = device.remove(duplicates)
Can you explain why set_nsc([1, 1, 1]) is called? And if so, is it necessary to do this step?
Ensure the lattice vectors are big enough for plotting.
Try and convince your-self that the lattice vectors are unimportant for tbtrans in this example.
HINT: what is the periodicity?
In [ ]:
device = device.add_vacuum(70, 0).add_vacuum(20, 1)
device = device.translate( device.center(what='cell') - device.center(what='xyz') )
device.write('device.xyz')
Since this system has 4 electrodes we need to tell tbtrans where the 4 electrodes are in the device. The following lines prints out the fdf-lines that are appropriate for each of the electrodes (RUN.fdf is already filled correctly):
In [ ]:
print('elec-Y-1: semi-inf -A2: {}'.format(1))
print('elec-Y-2: semi-inf +A2: end {}'.format(len(dev_y)))
print('elec-X-1: semi-inf -A1: {}'.format(len(dev_y) + 1))
print('elec-X-2: semi-inf +A1: end {}'.format(-1))
In [ ]:
H = sisl.Hamiltonian(device)
H.construct([R, hop])
H.write('DEVICE.nc')
In this example we have more than 1 transmission paths. Before you run the below code which plots all relevant transmissions ($T_{ij}$ for $j>i$), consider if there are any symmetries, and if so, determine how many different transmission spectra you should expect? Please plot the geometry using your favourite geometry viewer (molden, Jmol, ...). The answer is not so trivial.
In [ ]:
tbt = sisl.get_sile('siesta.TBT.nc')
Make easy function calls for plotting energy resolved quantites:
In [ ]:
E = tbt.E
Eplot = partial(plt.plot, E)
In [ ]:
# Make a shorthand version for the function (simplifies the below line)
T = tbt.transmission
t12, t13, t14, t23, t24, t34 = T(0, 1), T(0, 2), T(0, 3), T(1, 2), T(1, 3), T(2, 3)
Eplot(t12, label=r'$T_{12}$'); Eplot(t13, label=r'$T_{13}$'); Eplot(t14, label=r'$T_{14}$');
Eplot(t23, label=r'$T_{23}$'); Eplot(t24, label=r'$T_{24}$');
Eplot(t34, label=r'$T_{34}$');
plt.ylabel('Transmission'); plt.xlabel('Energy [eV]'); plt.ylim([0, None]); plt.legend();
RUN.fdf we have added the flag TBT.T.All which tells tbtrans to calculate all transmissions, i.e. between all $i\to j$ for all $i,j \in \{1,2,3,4\}$. This flag is by default False, why?
In [ ]:
# Insert plot of T12 and T21
In [ ]:
# Insert plot of T13 and T31
In [ ]:
# Insert plot of T14 and T41
In [ ]:
# Helper routines, this makes BDOS(...) == tbt.BDOS(..., norm='atom')
BDOS = partial(tbt.BDOS, norm='atom')
ADOS = partial(tbt.ADOS, norm='atom')
Bulk density of states:
In [ ]:
Eplot(..., label=r'$BDOS_1$');
Eplot(..., label=r'$BDOS_2$');
plt.ylabel('DOS [1/eV/N]'); plt.xlabel('Energy [eV]'); plt.ylim([0, None]); plt.legend();
Spectral density of states for all electrodes:
list) etc.
In [ ]:
Eplot(..., label=r'$ADOS_1$');
...
plt.ylabel('DOS [1/eV/N]'); plt.xlabel('Energy [eV]'); plt.ylim([0, None]); plt.legend();
matplotlib, here is the skeleton code for that, you should select an energy point and figure out how to extract the atom resolved DOS (you will need to look-up the documentation for the ADOS method to figure out which flag to use.
In [ ]:
Eidx = tbt.Eindex(...)
ADOS = [tbt.ADOS(i, ....) for i in range(4)]
f, axs = plt.subplots(2, 2, figsize=(10, 10))
a_xy = tbt.geometry.xyz[tbt.a_dev, :2]
for i in range(4):
A = ADOS[i]
A *= 100 / A.max() # normalize to maximum 100 (simply for plotting)
axs[i // 2][i % 2].scatter(a_xy[:, 0], a_xy[:, 1], A, c="bgrk"[i], alpha=.5);
plt.xlabel('x [Ang]'); plt.ylabel('y [Ang]'); plt.axis('equal');
In [ ]: