In [ ]:
from __future__ import print_function
import math
import tqdm
import sisl
import numpy as np
import matplotlib.pyplot as plt
from functools import partial
%matplotlib inline
In TBtrans and TranSiesta one is capable of performing real space transport calculations by using real space self-energies (see here).
Currently the real space self-energy calculation has to be performed in sisl since it is not implemented in TranSiesta.
The basic principle for calculating the real space self-energy is the Brillouin zone integral: \begin{equation} \mathbf G_{\mathcal R}(E) = \int_{\mathrm{BZ}}\mathbf G_\mathbf k \end{equation} In this example we will construct an STM tip probing a graphene flake.
We start by creating the graphene tight-binding model.
In [ ]:
graphene = sisl.geom.graphene(orthogonal=True)
In [ ]:
# Graphene tight-binding parameters
on, nn = 0, -2.7
H_minimal = sisl.Hamiltonian(graphene)
H_minimal.construct([[0.1, 1.44], [on, nn]])
Once the minimal graphene unit-cell (here orthogonal) is created we now turn to the calculation of the real space self-energy.
The construction of this object is somewhat complicated:
object: the Hamiltoniansemi_axes: which axes to use for the recursive self-energyk_axis: which axis to integrate in the Brillouin zoneunfold: how many times the object needs to be unfolded along each lattice vector, thus this is an integer vector of length 3
In [ ]:
RSSE = sisl.RealSpaceSE(H_minimal, 0, 1, (10, 10, 1))
Now we can create the real space self-energy.
However, in TBtrans (and TranSiesta) the electrode atomic indices must be in consecutive order.
This is a little troublesome since the natural order in a device would be an order according to $x$, $y$ or $z$. To create the correct order we extract the real space coupling matrix which is where the real space self-energy would like, the self-energy is calculated using:
\begin{equation}
\boldsymbol\Sigma^{\mathcal R} = E \mathbf S - \mathbf H - \Big[\int_{\mathrm{BZ}} \mathbf G\Big]^{-1}.
\end{equation}
Another way to calculate the self-energy would be to transfer the Green function from the infinite bulk into the region of interest:
\begin{equation}
\boldsymbol\Sigma^{\mathcal R} = \mathbf V_{\mathcal R\infty}\mathbf G_{\infty\setminus\mathcal R}\mathbf V_{\infty\mathcal R}.
\end{equation}
From the 2nd equation it is obvious that the self-energy only lives on the boundary that $\mathbf V_{\infty\mathcal R}$ couples to. Exactly this region is extracted using real_space_coupling as below.
In [ ]:
H_elec, elec_indices = RSSE.real_space_coupling(ret_indices=True)
H_elec.write('GRAPHENE.nc')
In [ ]:
H = RSSE.real_space_parent()
# Create the true device by re-arranging the atoms
indices = np.arange(len(H))
indices = np.delete(indices, elec_indices)
indices = np.concatenate([elec_indices, indices])
# Now re-arange matrix
H = H.sub(indices)
Lastly, we need to add the STM tip. Here we simply add a gold atom and manually add the hoppings.
In [ ]:
STM = sisl.Geometry([0, 0, 0], atom=sisl.Atom('Au', R=1.0001), sc=sisl.SuperCell([10, 10, 1], nsc=[1, 1, 3]))
H_STM = sisl.Hamiltonian(STM)
H_STM.construct([(0.1, 1.1), (0, -1)])
H_STM.write('STM.nc')
In [ ]:
mid_xyz = H.geometry.center()
idx = H.close(mid_xyz, R=1.33)[0]
H_device = H.add(H_STM, offset=H.geometry.xyz[idx] - H_STM.geometry.xyz[0] + [0, 0, 2])
na = len(H)
idx = H_device.close(na, R=(0.1, 2.25))[1][0]
H_device[na, idx] = -0.1
H_device[idx, na] = -0.1
In [ ]:
H_device.write('DEVICE.nc')
Before we can run calculations we need to create the real space self-energy for the graphene flake in sisl.
Since the algorithm is not implemented in TBtrans (nor TranSiesta) it needs to be done here.
This is somewhat complicated since the files requires a specific order. For ease this tutorial implements it for you.
In [ ]:
# A real space transport calculation ONLY needs the Gamma-point
gamma = sisl.MonkhorstPack(H_elec, [1] * 3)
# Energy contour
dE = 0.02
E = np.arange(-2, 2 + dE / 2, dE)
sisl.io.tableSile('contour.E', 'w').write_data(E, np.zeros(E.size) + dE)
In [ ]:
# Now create the file (should take around 3-4 minutes)
eta = 0.001 * 1j
with sisl.io.tbtgfSileTBtrans('GRAPHENE.TBTGF') as f:
f.write_header(gamma, E + eta)
for ispin, new_k, k, e in tqdm.tqdm(f):
if new_k:
f.write_hamiltonian(H_elec.Hk(format='array', dtype=np.complex128))
SeHSE = RSSE.self_energy(e + eta, bulk=True, coupling=True)
f.write_self_energy(SeHSE)
Calculate transport, density of state and bond-currents.
Before running the calculations you have to edit the input file to use the generated TBTGF file.
Additionally one have to edit the $E$ points used for transmission to use the contour.E file.
Please search the manual on how to edit the RUN.fdf file appropriately.
Plot the bond-currents and check their symmetry
TIME Redo the calculations using 3 electrodes (left/right/tip) using k-points. Converge transmission and then plot the bond-currents.
Do they look as the real space calculation?
In [ ]:
tbt = sisl.get_sile('siesta.TBT.nc')
In [ ]: