Initial commit of ts-tbt-sisl tutorial

This was the tutorial done in November 2016 in Barcelona

Signed-off-by: Nick Papior <nickpapior@gmail.com>

Author: zerothi

01/graphene_band_structure.py

@@ -0,0 +1,61 @@
+from __future__ import division
+import numpy as np
+
+def create_kpath(Nk):
+    # Gamma -> K -> M -> Gamma
+
+    # Gamma = 0   , 0
+    # K     = 2/3 , 1/3
+    # M     = 1/2 , 1/2
+
+    G2K = ( (2./3)**2 + (1./3)**2 ) **.5
+    K2M = ( (2./3 - 1./2)**2 + (1./3-1./2)**2 ) **.5
+    M2G = ( (1./2)**2 + (1./2)**2 ) **.5
+
+    Kdist = G2K + K2M + M2G
+
+    NG2K = int(Nk / Kdist * G2K)
+    NK2M = int(Nk / Kdist * K2M)
+    NM2G = int(Nk / Kdist * M2G)
+
+    def from_to(N, f, t):
+        full = np.empty([N, 3])
+        ls = np.linspace(0, 1, N, endpoint=False)
+        for i in range(3):
+            full[:, i] = f[i] + (t[i]-f[i]) * ls
+        return full
+
+    kG2K = from_to(NG2K, [0., 0., 0.], [2./3, 1./3, 0])
+    kK2M = from_to(NK2M, [2./3, 1./3, 0], [1./2, 1./2, 0.])
+    kM2G = from_to(NM2G, [1./2, 1./2, 0.], [0., 0., 0.])
+
+    xtick = [0, NG2K-1, NG2K + NK2M-1, NG2K + NK2M + NM2G-1]
+    label = ['G','K',         'M',                'G']
+
+    return [xtick, label], np.vstack((kG2K, kK2M, kM2G))
+
+def bandstructure(Nk, H):
+
+    ticks, k = create_kpath(Nk)
+    
+    eigs = np.empty([len(k),2],np.float64)
+    for ik, k in enumerate(k):
+        eigs[ik,:] = H.eigh(k=k)
+
+    import matplotlib.pyplot as plt
+
+    plt.plot(eigs[:,0])
+    plt.plot(eigs[:,1])
+    plt.gca().xaxis.set_ticks(ticks[0])
+    plt.gca().set_xticklabels(ticks[1])
+    ymin, ymax = plt.gca().get_ylim()
+    # Also plot x-major lines at the ticks
+    for tick in ticks[0]:
+        plt.plot([tick,tick], [ymin,ymax], 'k')
+    plt.show()
+        
+
+    
+if __name__ == "__main__":
+    print('This file is intended for use in Hancock_*.py files')
+    print('It is not intended to be runned separately.')

01/run.py

@@ -0,0 +1,111 @@
+#!/usr/bin/env python
+
+# This tutorial setup the necessary files for running
+# a simple graphene system example.
+
+# We will HEAVILY encourage to use Python3 compliant code
+# This line ensures the code will run using either Python2 or Python3
+from __future__ import print_function
+
+# First we require the use of the sisl Python module.
+import sisl
+
+# Instead of manually defining the graphene system we use
+# the build-in sisl capability of defining the graphene
+# structure.
+graphene = sisl.geom.graphene()
+
+# Print out some basic information about the geometry we have
+# just created.
+print("Graphene geometry information:")
+print(graphene)
+print()
+# It should print that the geometry consists of
+#   2 atoms (na)
+#   2 orbitals (no), one per atom
+#   1 specie (C), both atoms are the same atomic specie
+#   orbital range of 1.4342 AA (dR)
+#   3x3 supercell (nsc), this is determined from dR
+
+# Now as we have the graphene unit-cell, we can create the 
+# Hamiltonian for the system.
+H = sisl.Hamiltonian(graphene)
+print("Graphene initial Hamiltonian information:")
+print(H)
+print()
+
+# Currently this Hamiltonian is an empty entity, i.e.
+# no hopping elements has been assigned (all H[i,j] == 0).
+# For graphene, one may use these tight-binding parameters:
+#  on-site: 0. eV
+#  nearest-neighbour: -2.7 eV
+# In the following loop we setup the Hamiltonian with the above
+# tight-binding parameters:
+for ia, io in H:
+	# This loop construct loops over all atoms and the orbitals
+    # corresponding to the atom.
+	# In this case the geometry has one orbital per atom, hence
+    #   ia == io
+    # in all cases.
+    
+    # In order to figure out which atoms atom `ia` is connected
+    # to, we must find those atoms
+    # To do this we access the geometry attached to the 
+    # Hamiltonian (H.geom)
+    # and use a function called `close` which returns ALL 
+    # atoms within certain ranges of a given point or atom
+    idx = H.geom.close(ia, dR = (0.1, 1.43))
+    # the argument dR has two entries:
+    #   0.1 and 1.43
+    # Each value represents a radii of a sphere.
+    # The `close` function will then return
+    # a list of equal length of the dR argument (i.e. a list with
+    # two values).
+    # Then idx[0] is the first element and this contains a list
+    # of all atoms within a sphere of 0.1 AA of atom `ia`.
+    # This should obviously only contain the atom it-self.
+    # The second element, idx[1] contains all atoms within a sphere
+    # with radius of 1.43 AA, but not including those within 0.1 AA.
+    # In this case this is then all atoms that are the nearest neighbour
+    # atoms
+   
+    # Now we know the on-site atoms (idx[0]) and the nearest neighbour
+    # atoms (idx[1]), all we need to do is set the Hamiltonian
+    # elements:
+
+    # on-site (0. eV)
+    H[io, idx[0]] = 0.
+   
+    # nearest-neighbour (-2.7 eV)
+    H[io, idx[1]] = -2.7
+    
+# At this point we have created all Hamiltonian elements for all orbitals
+# in the graphene geometry.
+# Let's try and print the information of the Hamiltonian object:
+print("Hamiltonian after setting up the parameters:")
+print(H)
+print()
+# It now prints that there are 8 non-zero elements.
+# Please convince yourself that this is the correct number of 
+# Hamiltonian elements (regard the on-site value of 0. eV as a non-zero
+# value). HINT: count the number of on-site and nearest-neighbour terms)
+
+# At this point we have all the ingredients for manipulating the electronic
+# structure of graphene in a tight-binding model with nearest neighbour 
+# interactions.
+
+# We may, for instance, calculate the eigen-spectrum at the Gamma-point:
+print("Eigenvalues at Gamma:")
+print(H.eigh())
+print("Eigenvalues at K:")
+print(H.eigh(k=[1./3,2./3,0]))
+print()
+
+# Additionally we may calculate the band-structure of graphene
+# We want to plot the band-structure of graphene
+
+from graphene_band_structure import bandstructure
+
+# Now lets plot the band-structure (requires matplotlib)
+bandstructure(301, H)
+

02/RUN.fdf

@@ -0,0 +1,46 @@
+# Define the device Hamiltonian
+TBT.HS DEVICE.nc
+
+# Graphene requires a high k-point sampling
+# Try and play with this number to find
+# a value which gives a smooth
+# integration.
+# Please note which unit-cell direction
+# that the below array defines
+TBT.k [50 1 1]
+
+# These are the definitions of the electrodes
+# They may seem very verbose, but this
+# gives the user a very high degree of
+# flexibility. This will be apparent
+# in later tutorials
+%block TBT.Elec.Left
+  # Define the electronic structure
+  # of the electrode
+  HS ELEC.nc
+
+  # The semi-infinite direction was
+  # along the second lattice vector (A2)
+  # The negative sign is to specify 
+  # the direction of the semi-infinite
+  # bulk part
+  semi-inf-direction -A2
+
+  # Denote the first atom of the electrode.
+  # This index corresponds to the atomic
+  # index in the TBT.HS file
+  electrode-position 1
+%endblock TBT.Elec.Left
+
+# Also define the equivalent quantities 
+# for the right electrode
+%block TBT.Elec.Right
+  HS ELEC.nc
+  # The system is <- Left | device | Right ->
+  # hence the semi-infinite direction of the Right
+  # electrode has +
+  semi-inf-direction +A2
+  # Note the use of "end" which 
+  # refers to the last atom of the electrode
+  electrode-position end -1
+%endblock TBT.Elec.Right

02/run.py

@@ -0,0 +1,50 @@
+#!/usr/bin/env python
+
+import sisl
+
+# Instead of manually defining the graphene system we use
+# the build-in sisl capability of defining the graphene
+# structure.
+# Note that this graphene unit-cell is an orthogonal unit-cell
+# (figure out how many atoms this consists off)
+graphene = sisl.geom.graphene(orthogonal=True)
+
+# Now as we have the graphene unit-cell, we can create the 
+# Hamiltonian for the system.
+H = sisl.Hamiltonian(graphene)
+
+# Instead of using the for-loop provided in 01, we
+# may use a simpler function that *ONLY* works
+# for 1-orbital per atom models.
+H.construct([0.1, 1.43], [0., -2.7])
+    
+# At this point we have created all Hamiltonian elements for all orbitals
+# in the graphene geometry.
+# Please try and figure out how many non-zero elements there should
+# be in the Hamiltonian object.
+# HINT: count the number of on-site and nearest-neighbour terms
+
+# Now we have the Hamiltonian, now we need to save it to be able 
+# to conduct a tbtrans calculation:
+# As this is the minimal unit-cell we call this the electrode
+# and we will in the following create a larger "DEVICE" region
+# Hamiltonian.
+H.write('ELEC.nc')
+
+# Now we need to create the device region Hamiltonian.
+# The tbtrans method relies on including the electrodes 
+# in the device region. Hence the smallest device region 
+# must be 3 times the electrode size.
+# Please try and convince your-self that this is the case.
+# HINT: consider the interaction range of the atoms
+
+# We tile the geometry along the 2nd lattice vector
+# (Python is 0-based indexing)
+device = graphene.tile(3, axis=1)
+Hdev = sisl.Hamiltonian(device)
+Hdev.construct([0.1, 1.43], [0, -2.7])
+# We will create both a xyz file (for plotting in molden, etc.)
+# and we will create the Hamilton file format for reading in tbtrans
+Hdev.geom.write('device.xyz')
+Hdev.write('DEVICE.nc')
+

03/RUN.fdf

@@ -0,0 +1,30 @@
+# Define the device Hamiltonian
+TBT.HS DEVICE.nc
+
+# Graphene requires a high k-point sampling
+# Try and play with this number to find
+# a value which gives a smooth
+# integration.
+TBT.k [50 1 1]
+
+# These are the definitions of the electrodes
+# They may seem very verbose, but this
+# gives the user a very high degree of
+# flexibility. This will be apparent
+# in later tutorials
+%block TBT.Elec.Left
+  HS ELEC.nc
+  # The semi-infinite direction was
+  # along the second lattice vector (A2)
+  # The negative sign is to specify 
+  # the direction of the semi-infinite
+  # bulk part
+  semi-inf-direction -A2
+  electrode-position 1
+%endblock TBT.Elec.Left
+
+%block TBT.Elec.Right
+  HS ELEC.nc
+  semi-inf-direction +A2
+  electrode-position end -1
+%endblock TBT.Elec.Right

03/run.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python
+
+import sisl
+
+# Instead of manually defining the graphene system we use
+# the build-in sisl capability of defining the graphene
+# structure.
+# (figure out how many atoms this consists off)
+graphene = sisl.geom.graphene().tile(2, axis=0)
+
+# Now as we have the graphene unit-cell, we can create the 
+# Hamiltonian for the system.
+H = sisl.Hamiltonian(graphene)
+
+# Instead of using the for-loop provided in 01, we
+# may use a simpler function that *ONLY* works
+# for 1-orbital per atom models.
+H.construct([0.1, 1.43], [0., -2.7])
+    
+# At this point we have created all Hamiltonian elements for all orbitals
+# in the graphene geometry.
+# Please try and figure out how many non-zero elements there should
+# be in the Hamiltonian object.
+# HINT: count the number of on-site and nearest-neighbour terms
+
+# Now we have the Hamiltonian, now we need to save it to be able 
+# to conduct a tbtrans calculation:
+# As this is the minimal unit-cell we call this the electrode
+# and we will in the following create a larger "DEVICE" region
+# Hamiltonian.
+H.write('ELEC.nc')
+
+# Now we need to create the device region Hamiltonian.
+# The tbtrans method relies on including the electrodes 
+# in the device region. Hence the smallest device region 
+# must be a multiple of 3 of the electrode size.
+# Please try and convince your-self that this is the case.
+# HINT: consider the interaction range of the atoms
+
+# We tile the geometry along the 2nd lattice vector
+# (Python is 0-based indexing)
+device = graphene.tile(3, axis=1)
+Hdev = sisl.Hamiltonian(device)
+Hdev.construct([0.1, 1.43], [0, -2.7])
+# We will create both a xyz file (for plotting in molden, etc.)
+# and we will create the Hamilton file format for reading in tbtrans
+Hdev.geom.write('device.xyz')
+Hdev.write('DEVICE.nc')
+

04/RUN.fdf

@@ -0,0 +1,32 @@
+TBT.HS DEVICE.nc
+
+# This is also a transverse periodic system and
+# again the k-mesh is important for converged results.
+TBT.k [5 1 1]
+
+# Define which physical quantites to calculate.
+# Each of these flags will decide what to calculate
+# in the calculation.
+# If you look these flags up in the manual you will
+# find the physical quantities they correspond to.
+TBT.DOS.Gf true
+TBT.DOS.Elecs true
+TBT.DOS.A true
+TBT.DOS.A.All true
+TBT.T.All true
+TBT.T.Bulk true
+TBT.T.Eig 5
+TBT.Current.Orb true
+TBT.Symmetry.TimeReversal false
+
+
+%block TBT.Elec.Left
+  HS ELEC.nc
+  semi-inf-direction -A2
+  electrode-position 1
+%endblock TBT.Elec.Left
+%block TBT.Elec.Right
+  HS ELEC.nc
+  semi-inf-direction +A2
+  electrode-position end -1
+%endblock TBT.Elec.Right