fixed examples

examples/bi_bulk.py

@@ -1,5 +1,7 @@
 """
-This script computes band structure of the crystalline bismuth
+This example script computes band structure of the crystalline bismuth.
+It uses the third-nearest neighbor approximation with a step-wise distance distance
+associating distances with sets of TB parameters.
 """
 import numpy as np
 import matplotlib.pyplot as plt
@@ -13,7 +15,6 @@ def radial_dep(coords):
     """
         Step-wise radial dependence function
     """
-
     norm_of_coords = np.linalg.norm(coords)
     if norm_of_coords < 3.3:
         return 1

examples/nanostrip.py

@@ -1,6 +1,6 @@
 """
-This script computes DOS and transmission function for the four-atom-width nanostrip
-using the recursive Green's function algorithm
+This example script computes DOS and transmission function for the four-atom-width nanostrip
+using the recursive Green's function algorithm.
 """
 import numpy as np
 import nanonet.tb as tb
@@ -9,7 +9,6 @@
 
 
 def main(energy):
-
     # define the basis set - one s-type orbital
     orb = tb.Orbitals('A')
     orb.add_orbital('s', energy=-1.0)

examples/si_bulk.py

@@ -1,3 +1,6 @@
+"""
+This example script computes band structure of the bulk silicon.
+"""
 import numpy as np
 import matplotlib.pyplot as plt
 from nanonet.tb import Hamiltonian
@@ -6,19 +9,14 @@
 
 
 # Crystal structure parameters
-
 a = 5.50
 
-
 # Primitive cell
-
 primitive_cell = a * np.array([[0.0, 0.5, 0.5],
                                [0.5, 0.0, 0.5],
                                [0.5, 0.5, 0.0]])
 
-
 # High symmetry points
-
 SPECIAL_K_POINTS_SI = {
     'GAMMA': [0, 0, 0],
     'X': [0, 2 * np.pi / a, 0],
@@ -30,24 +28,34 @@
 
 
 def main():
+    # specify atomic coordinates in xyz format
+    path_to_xyz_file = """2
+                          Bulk Si cell
+                          Si1       0.000    0.000    0.000
+                          Si2       1.375    1.375    1.375"""
 
-    path_to_xyz_file = 'input_samples/bulk_silicon.xyz'
+    # specify basis set
     Orbitals.orbital_sets = {'Si': 'SiliconSP3D5S'}
 
-    sym_points = ['L', 'GAMMA', 'X']
-    num_points = [20, 20]
-
+    # create a Hamiltonian object storing the Hamiltonian matrices
     h = Hamiltonian(xyz=path_to_xyz_file)
     h.initialize()
+
+    # set periodic boundary conditions
     h.set_periodic_bc(primitive_cell)
 
+    # define wave vector coordinates
+    sym_points = ['L', 'GAMMA', 'X']
+    num_points = [20, 20]
     k_points = get_k_coords(sym_points, num_points, SPECIAL_K_POINTS_SI)
 
+    # compute band structure
     band_structure = []
     for ii, item in enumerate(k_points):
         eigenvalues, _ = h.diagonalize_periodic_bc(k_points[ii])
         band_structure.append(eigenvalues)
 
+    # visualize
     band_structure = np.array(band_structure)
     ax = plt.axes()
     ax.plot(band_structure)

examples/si_nw_band_structure.py

@@ -15,7 +15,7 @@ def main():
     # specify atomic coordinates file in xyz format
     path = "input_samples/SiNW2.xyz"
 
-    # define Hamiltonian object in the sparse format
+    # create a Hamiltonian object in the sparse format
     hamiltonian = HamiltonianSp(xyz=path, nn_distance=2.4,)
     hamiltonian.initialize()
 

examples/si_nw_transmission.py

@@ -1,88 +1,90 @@
+"""
+This example script computes DOS and transmission function for the silicon nanowire
+using the recursive Green's function algorithm.
+"""
 import logging
-import time
 import numpy as np
 import nanonet.negf as negf
 import matplotlib.pyplot as plt
-from nanonet.tb import Hamiltonian, HamiltonianSp
+from nanonet.tb import Hamiltonian
 from nanonet.tb import Orbitals
 from nanonet.tb.sorting_algorithms import sort_projection
 
 
-def radial_dep(coords):
-    norm_of_coords = np.linalg.norm(coords)
-    print(norm_of_coords)
-    if norm_of_coords < 3.3:
-        return 1
-    elif 3.7 > norm_of_coords > 3.3:
-        return 2
-    elif 5.0 > norm_of_coords > 3.7:
-        return 3
-    else:
-        return 100
-
-
 def main():
-    # define orbitals sets
+    # use a predefined basis sets
     Orbitals.orbital_sets = {'Si': 'SiliconSP3D5S', 'H': 'HydrogenS'}
 
+    # specify atomic coordinates file stored in xyz format
     path = "input_samples/SiNW2.xyz"
 
+    # define leads indices
     right_lead = [0, 33, 35, 17, 16, 51, 25,  9, 53, 68,  1,  8, 24]
     left_lead = [40, 66, 58, 47, 48, 71, 72, 73, 74, 65]
 
-    start = time.time()
-    hamiltonian = Hamiltonian(xyz=path, nn_distance=2.4, sort_func=sort_projection, left_lead=left_lead, right_lead=right_lead)
+    # create a Hamiltonian object storing the Hamiltonian matrices
+    hamiltonian = Hamiltonian(xyz=path, nn_distance=2.4,
+                              sort_func=sort_projection,
+                              left_lead=left_lead, right_lead=right_lead)
     hamiltonian.initialize()
 
+    # set periodic boundary conditions
     a_si = 5.50
     primitive_cell = [[0, 0, a_si]]
     hamiltonian.set_periodic_bc(primitive_cell)
 
-    forming_time = time.time() - start
-    print(forming_time)
-
+    # get Hamiltonian matrices
     hl, h0, hr = hamiltonian.get_hamiltonians()
+    # get Hamiltonian matrices in the block-diagonal form
     hl1, h01, hr1, subblocks = hamiltonian.get_hamiltonians_block_tridiagonal()
 
-    energy = np.linspace(2.1, 2.5, 50)
-    tr = np.zeros(energy.shape)
-    dos = np.zeros(energy.shape)
+    # specify energy array
+    energy = np.linspace(2.1, 2.3, 30)
 
+    # specify dephasing constant
     damp = 0.0005j
 
+    # initialize output arrays by zeros
+    tr = np.zeros(energy.shape)
+    dos = np.zeros(energy.shape)
+
+    # energy loop
     for j, E in enumerate(energy):
+
         logging.info("{} Energy: {} eV".format(j, E))
-        start = time.time()
+
+        # compute self-energies describing boundary conditions at the leads contacts
         R, L = negf.surface_greens_function(E, hl, h0, hr, iterate=True, damp=damp)
-        sgf_time = time.time() - start
-        print(sgf_time)
+
         s01, s02 = h01[0].shape
         s11, s12 = h01[-1].shape
 
+        # apply self-energies
         h01[0] = h01[0] + L[:s01, :s02]
         h01[-1] = h01[-1] + R[-s11:, -s12:]
 
-        start = time.time()
+        # compute Green's functions using the recursive Green's function algorithm
         # g_trans, grd, grl, gru, gr_left = negf.recursive_gf(E, [hl, hl], [h0+L, h0, h0+R], [hr, hr], damp=damp)
         g_trans, grd, grl, gru, gr_left = negf.recursive_gf(E, hl1, h01, hr1, damp=damp)
-        rgf_time = time.time() - start
-        print(rgf_time)
 
+        # detach self-energies
         h01[0] = h01[0] - L[:s01, :s02]
         h01[-1] = h01[-1] - R[-s11:, -s12:]
 
-        num_periods = 3
+        # number of subblocks
+        num_periods = len(grd)
 
+        # compute DOS
         for jj in range(num_periods):
             dos[j] = dos[j] + np.real(np.trace(1j * (grd[jj] - grd[jj].conj().T))) / num_periods
 
-        # dos[j] = -2 * np.trace(np.imag(grd[0]))
-
         gamma_l = 1j * (L[:s01, :s02] - L[:s01, :s02].conj().T)
         gamma_r = 1j * (R[-s11:, -s12:] - R[-s11:, -s12:].conj().T)
 
+        # compute transmission spectrum
         tr[j] = np.real(np.trace(gamma_l.dot(g_trans).dot(gamma_r).dot(g_trans.conj().T)))
 
+    # visualize
     plt.plot(energy, dos)
     plt.show()
 

nanonet/negf/greens_functions.py

@@ -7,6 +7,7 @@
 import os.path
 import numpy as np
 import scipy.linalg as linalg
+from scipy.sparse import csr_matrix
 
 
 def surface_greens_function_poles(h_list):