In [ ]:
from __future__ import print_function
import math
import sisl
import numpy as np
import matplotlib.pyplot as plt
from functools import partial
%matplotlib inline
One of the key ideas behind sisl was the interaction of a DFT Hamiltonian and the user.
In this example we will highlight a unique implementation in TBtrans which enables any kind of user intervention.
The idea is a transformation of the Green function calculation from:
\begin{equation}
\mathbf G^{-1}(E) = \mathbf S (E + i\eta) - \mathbf H - \sum_i\boldsymbol\Sigma_i
\end{equation}
to
\begin{equation}
\mathbf G^{-1}(E) = \mathbf S (E + i\eta) - \mathbf H - \delta\mathbf H - \sum_i\boldsymbol\Sigma_i - \delta\boldsymbol \Sigma
\end{equation}
where $\delta\mathbf H$ and $\delta\boldsymbol\Sigma$ can be of any type, i.e. complex and/or real.
The only (important!) difference between $\delta\mathbf H$ and $\delta\boldsymbol \Sigma$ is that the former enters the calculation of bond-currents, while the latter does not.
Since TBtrans by it-self does not allow complex Hamiltonians the above is a way to leviate this restriction. One feature this may be used for is by applying magnetic fields.
In the following we will use Peierls substitution on a square tight-binding model: \begin{equation} \mathbf H(\Phi) = \mathbf H(0)e^{i \Phi \delta x^- \cdot \delta y^+ / 2}, \end{equation} where $\Phi$ is the magnetic flux (in proper units), $\delta x^- = x_j - x_i$ and $\delta y^+=y_j + y_i$, with $i$ and $j$ being atomic indices. If you are interested in this substitution, please search litterature.
First create a square lattice and define the on-site and nearest neighbour couplings
In [ ]:
square = sisl.Geometry([0,0,0], sisl.Atom(1, R=1.0), sc=sisl.SuperCell(1, nsc=[3, 3, 1]))
on, nn = 4, -1
In [ ]:
H_minimal = sisl.Hamiltonian(square)
H_minimal.construct([[0.1, 1.1], [on, nn]])
H_elec = H_minimal.tile(100, 1).tile(2, 0)
H_elec.set_nsc([3, 1, 1])
H_elec.write('ELEC.nc')
In [ ]:
H = H_elec.tile(50, 0)
# Make a constriction
geom = H.geom.translate( -H.geom.center(what='xyz') )
# This constriction is based on an example in the kwant project (called qhe). We however make a slight modification.
# Remove some atoms, this will create a constriction of 100 - 40 * 2 = 20 Ang with a Gaussian edge profile
remove = (np.abs(geom.xyz[:, 1]) > 50 - 37.5 * np.exp( -(geom.xyz[:, 0] / 12) **2 )).nonzero()[0]
# To reduce computations we find the atoms in the constriction such that we can
# limit the calculation region.
device = (np.abs(geom.remove(remove).xyz[:, 0]) < .6).nonzero()[0]
geom.remove(remove).write('test.xyz')
# Pretty print a range of atoms that is the smallest device region
print(sisl.utils.list2str(device))
H = H.remove(remove)
H.write('DEVICE.nc')
The above printed list of atoms should be inserted in the RUN.fdf in the TBT.Atoms.Device. This is important when one is only interested in the transmission, and does not care about the density of states. You are always encouraged to select the minimal device region to 1) speed up computations and 2) drastically reduce memory requirements.
In this regard you should read in the TBtrans manual about the additional flag TBT.Atoms.Device.Connect (you may try, as an additional exercise to set this flag to true and check the difference from the previous calculation).
Now we have $\mathbf H(0)$ with no phases due to magnetic fields. As the magnetic field is changing the Hamiltonian, and thus enters the bond-current calculations, we have to use the $\delta\mathbf H$ term (and not $\delta\boldsymbol\Sigma$).
The first thing we need to calculate is $\delta x^- \cdot\delta y^+$. Since we already have the Hamiltonian we can utilize the connections by looping the coupling elements (in a sparse matrix/graph this is called edges). This is much cheaper than trying to figure out which atoms are neighbouring.
In [ ]:
device = H.geom
xy = sisl.Hamiltonian(device)
for ia in device:
# Get all connecting elements (excluding it-self)
edges = H.edges(ia)
# Calculate the vector between edges and ia:
# xyz[edges, :] - xyz[ia, :]
Rij = device.Rij(ia, edges)
# Now calculate the product:
# (xj - xi) * (yj + yi)
# Notice that we correct to +yi by adding it twice
xy[ia, edges] = Rij[:, 0] * (Rij[:, 1] + 2 * device.xyz[ia, 1])
xy[ia, ia] = 0.
xy.finalize() # this is only because it will speed up the following calculations
Now we have the coupling dependent phase factor, $\delta x^-\cdot\delta y^+$, and all we need to calculate is $\delta\mathbf H$ that transforms $\mathbf H(0) \to \mathbf H(\Phi)$.
This is done easily as a Hamiltonian allows basic element wise operations, i.e. +, -, *, / and ** (the power function).
Your task is to insert the correct mathematical equation below, such that dH contains $\delta \mathbf H$ for $\mathbf H(\Phi) = \mathbf H(0) + \delta\mathbf H$.
To help you I have inserted the exponential function. To finalize the equation, you need three terms: nn, xy and rec_phi.
In [ ]:
rec_phis = np.arange(1, 51, 4)
print('Calulating (of {}):'.format(len(rec_phis)), end='')
for i, rec_phi in enumerate(rec_phis):
print(' {}'.format(i+1), end=',')
# Calculate H(Phi)
# TODO: insert the correct mathematical equation below to calculate the correct phase
#dH = math.e ** (0.5j ... ) ...
with sisl.get_sile('M_{}.dH.nc'.format(rec_phi), mode='w') as fh:
fh.write_delta(dH)
Calculate all physical quantities for all different applied magnetic fields.
Before running the calculations, search the manual on how to save the self-energies (HINT out-of-core). By default, TBtrans calculates the self-energies as they are needed. However, if one has the same electrodes, same $k$-grid and same $E$-points for several different runs (as in this case) one can with benefit calculate the self-energies once, and then reuse them in subsequent calculations.
To ease the calculation of all magnetic fields
To help you a script run.sh is located in this directory. Carefully read it to infer which option specifies the $\delta \mathbf H$ term.
Since this example has 14 different setups, each with 51 energy points, it will take some time. Around 30 seconds for the first (includes self-energy calculation), and around 10 seconds for all subsequent setups. So be patient. :)
In [ ]:
# Create short-hand function
gs = sisl.get_sile
# No magnetic field
tbt0 = gs('siesta.TBT.nc')
# All magnetic fields in increasing order
tbts = [gs('M_{}/siesta.TBT.nc'.format(rec_phi)) for rec_phi in rec_phis]
In [ ]:
# Do trick with easy plotting utility
E = tbt0.E[:]
Eplot = partial(plt.plot, E)
for rec_phi, tbt in zip(rec_phis, tbts):
Eplot(tbt.transmission(), '--', label='{}'.format(rec_phi));
Eplot(tbt0.transmission(), 'k');
plt.xlim([E.min(), E.max()]); plt.ylim([0, None]);
plt.xlabel('Energy [eV]'); plt.ylabel('Transmisson'); plt.legend(bbox_to_anchor=(1, 1), loc=2);
In [ ]:
T = np.stack([tbt.transmission() for tbt in tbts])
T[T < 1e-9] = 0 # small numbers are difficult to interpolate (they tend to blow up), so remove them
plt.contourf(rec_phis, E, T.T, 20); plt.colorbar(label='Transmission');
plt.ylabel('Energy [eV]'); plt.xlabel(r'$1/\Phi$');
In [4-7] and adapt), you may decide the Gaussian profile and change if you want.
In [ ]:
In [ ]:
TBT.dH vs. TBT.dSE). Please look in the manual and figure out what the difference is between the two methods?deltancSileTBtrans.
In [ ]: