In [ ]:
from __future__ import print_function
import sisl
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In this example we will begin to cover some of the extraction utilities in sisl allowing one to really go in-depth on analysis of calculations.
We will begin by creating a large graphene flake, then subsequently a hole will be created by removing a circular shape.
Subsequent to the calculations you are encouraged to explore the sisl toolbox for ways to extract important information regarding your system.
In [ ]:
graphene = sisl.geom.graphene(orthogonal=True)
In [ ]:
elec = graphene.tile(25, axis=0)
H = sisl.Hamiltonian(elec)
H.construct(([0.1, 1.43], [0., -2.7]))
H.write('ELEC.nc')
In [ ]:
device = elec.tile(15, axis=1)
device = device.remove(device.close(device.center(what='cell'), R=10.))
We will also (for physical reasons) remove all dangling bonds, and secondly we will create a list of atomic indices which corresponds to the atoms at the edge of the hole.
In [ ]:
dangling = [ia for ia in device.close(device.center(what='cell'), R=14.)
if len(device.close(ia, R=1.43)) < 3]
device = device.remove(dangling)
edge = []
for ia in device.close(device.center(what='cell'), R=14.):
if len(device.close(ia, R=1.43)) < 4:
edge.append(ia + 1) # + 1 to get fortran indices
# Pretty-print the list of atoms
print(sisl.utils.list2str(edge))
Hdev = sisl.Hamiltonian(device)
Hdev.construct(([0.1, 1.43], [0, -2.7]))
Hdev.geom.write('device.xyz')
Hdev.write('DEVICE.nc')
Please carefully go through the RUN.fdf file before running TBtrans, determine what each flag means and what it tells TBtrans to output.
Now run TBtrans:
tbtrans RUN.fdf
In addition to the previous example, many more files are now being created (for all files the siesta.TBT.AV* files are the $k$-averaged equivalents. You should, while reading the below list, also be able to specify which of the fdf-flags that are responsible for the creation of which files.
siesta.TBT.ADOS_*siesta.TBT.BDOS_*siesta.TBT.BTRANS_*siesta.TBT.DOSsiesta.TBT.TEIG_<1>-<2><1> and emitted in <2>.This exercise mainly shows a variety of methods useful for extracting data from the *.TBT.nc files in a simple and consistent manner. You are encouraged to play with the routines, because the next example forces you to utilize them.
In [ ]:
tbt = sisl.get_sile('siesta.TBT.nc')
As this system is not a pristine periodic system we have a variety of options available for analysis.
First and foremost we plot the transmission:
In [ ]:
plt.plot(tbt.E, tbt.transmission(),label='Av');
plt.plot(tbt.E, tbt.transmission(kavg=tbt.kindex([0,0,0])), label=r'$\Gamma$'); plt.legend()
plt.ylabel('Transmission'); plt.xlabel('Energy [eV]'); plt.ylim([0, None]);
The contained data in the *.TBT.nc file is very much dependent on the flags in the fdf-file. To ease the overview of the available output and what is contained in the file one can execute the following block to see the content of the file.
Check whether the bulk transmission is present in the output file and if so, add it to the plot above to compare the bulk transmission with the device transmission.
There are two electrodes, hence two bulk transmissions. Is there a difference between the two? If so, why, if not, why not?
In [ ]:
print(tbt.info())
We may also plot the Green function density of states as well as the spectral density of states:
In [ ]:
DOS_all = [tbt.DOS(), tbt.ADOS(0), tbt.ADOS(1)]
plt.plot(tbt.E, DOS_all[0], label='G');
plt.plot(tbt.E, DOS_all[1], label=r'$A_L$');
plt.plot(tbt.E, DOS_all[2], label=r'$A_R$');
plt.ylim([0, None]); plt.ylabel('Total DOS [1/eV]'); plt.xlabel('Energy [eV]'); plt.legend();
Can you from the above three quantities determine whether there are any localized states in the examined system?
HINT: What is the sum of the spectral density of states ($\mathbf A$) compared to the Green function ($\mathbf G$) density of states?
The total DOS is a measure of the DOS spread out in the entire atomic region. However, TBtrans calculates, and stores all DOS related quantities in orbital resolution. I.e. we are able to post-process the DOS and examine the atom (orbital) resolved DOS.
To do this the .DOS and .ADOS routines have two important keywords, 1) atom and 2) orbital which may be used to extract a subset of the DOS quantities. For details on extracting these subset quantities please read the documentation by executing the following line:
help(tbt.DOS)
The following code will extract the DOS only on the atoms in the hole edge.
In [ ]:
DOS_edge = [tbt.DOS(atom=edge), tbt.ADOS(0, atom=edge), tbt.ADOS(1, atom=edge)]
plt.plot(tbt.E, DOS_edge[0], label='G');
plt.plot(tbt.E, DOS_edge[1], label=r'$A_L$');
plt.plot(tbt.E, DOS_edge[2], label=r'$A_R$');
plt.ylim([0, None]); plt.ylabel('DOS on edge atoms [1/eV]'); plt.xlabel('Energy [eV]'); plt.legend();
Comparing the two previous figures leaves us with little knowlegde of the DOS ratio. I.e. both plots show the total DOS and they are summed over a different number of atoms (or orbitals if you wish). Instead of showing the total DOS we can normalize the DOS by the number of atoms; $1/N_a$ where $N_a$ is the number of atoms in the selected DOS region. With this normalization we can compare the average DOS on all atoms with the average DOS on only edge atoms.
The tbtncSileTBtrans can readily do this for you.
Read about the norm keyword in .DOS, and also look at the documentation for the .norm function:
help(tbt.DOS)
help(tbt.norm)
Now create a plot with DOS normalized per atom by editing the below lines, feel free to add the remaining DOS plots to have them all:
In [ ]:
N_all = tbt.norm(<change here>)
N_edge = tbt.norm(<change here>)
plt.plot(tbt.E, DOS_all[0] / N_all, label=r'$G$');
plt.plot(tbt.E, DOS_edge[0] / N_edge, label=r'$G_E$');
plt.ylim([0, None]); plt.ylabel('DOS [1/eV/atom]'); plt.xlabel('Energy [eV]'); plt.legend();
We can further analyze the DOS evolution for atoms farther and farther from the hole.
The following DOS analysis will extract DOS (from $\mathbf G$) for the edge atoms, then for the nearest neighbours to the edge atoms, and for the next-nearest neighbours to the edge atoms.
Try and extend the plot to contain the DOS of the next-next-nearest neighbours to the edge atoms.
In [ ]:
# Get nearest neighbours to the edge atoms
n_edge = Hdev.edges(edge)
# Get next-nearest neighbours to the edge atoms
nn_edge = Hdev.edges(n_edge, exclude=np.hstack((edge, n_edge)))
# Try and create the next-next-nearest neighbours to the edge atoms and add it to the plot
plt.plot(tbt.E, tbt.DOS(atom=edge, norm='atom'), label='edge: G');
plt.plot(tbt.E, tbt.DOS(atom=n_edge, norm='atom'), label='n-edge: G');
plt.plot(tbt.E, tbt.DOS(atom=nn_edge, norm='atom'), label='nn-edge: G');
plt.ylim([0, None]); plt.ylabel('DOS [1/eV/atom]'); plt.xlabel('Energy [eV]'); plt.legend();