graphene = sisl.geom.graphene(orthogonal=True)

elec = graphene.tile(25, axis=0)
H = sisl.Hamiltonian(elec)
H.construct(([0.1, 1.43], [0., -2.7]))
H.write('ELEC.nc')

device = elec.tile(15, axis=1)
device = device.remove(
    device.close(
        device.center(what='cell'), R=10.)
)

"""
Explanation: 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.
End of explanation
"""


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)
edge = np.array(edge)

# Pretty-print the list of atoms (for use in sdata)
# Note we add 1 to get fortran indices
print(sisl.utils.list2str(edge + 1))

Hdev = sisl.Hamiltonian(device)
Hdev.construct(([0.1, 1.43], [0, -2.7]))
Hdev.geometry.write('device.xyz')
Hdev.write('DEVICE.nc')

"""
Explanation: 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.
End of explanation
"""


tbt = sisl.get_sile('siesta.TBT.nc')

"""
Explanation: Exercises
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_*
   The $k$-resolved spectral density of states injected from the named electrode
siesta.TBT.BDOS_*
   The $k$-resolved bulk Green function density of states for the named electrode
siesta.TBT.BTRANS_*
   The $k$-resolved bulk transmission through the named electrode
siesta.TBT.DOS
   Green function density of states
siesta.TBT.TEIG_&lt;1&gt;-&lt;2&gt;
   Transmission eigenvalues for electrons injected from &lt;1&gt; and emitted in &lt;2&gt;.

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.
End of explanation
"""


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]);

"""
Explanation: 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:
End of explanation
"""


print(tbt.info())

"""
Explanation: 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?
End of explanation
"""


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();

"""
Explanation: Density of states
We may also plot the Green function density of states as well as the spectral density of states from each electrode:
End of explanation
"""


DOS_edge = [tbt.DOS(atoms=edge), tbt.ADOS(0, atoms=edge), tbt.ADOS(1, atoms=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();

"""
Explanation: 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_i$) compared to the Green function ($\mathbf G$) density of states?
Examining DOS on individual atoms
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 and/or 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.
End of explanation
"""


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();

"""
Explanation: 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:
End of explanation
"""


# Get nearest neighbours to the edge atoms
n_edge = Hdev.edges(edge, exclude=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(atoms=edge, norm='atom'), label='edge: G');
plt.plot(tbt.E, tbt.DOS(atoms=n_edge, norm='atom'), label='n-edge: G');
plt.plot(tbt.E, tbt.DOS(atoms=nn_edge, norm='atom'), label='nn-edge: G');
plt.ylim([0, None]); plt.ylabel('DOS [1/eV/atom]'); plt.xlabel('Energy [eV]'); plt.legend();

"""
Explanation: DOS depending on the distance from the hole
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.
End of explanation
"""


