refactored greens_functions.py

.gitattributes

@@ -0,0 +1 @@
+jupyter-notebooks/* linguist-documentation

examples/si_nw_transmission_blocks.py

@@ -0,0 +1,93 @@
+"""
+This example script computes DOS and transmission function for the silicon nanowire
+using the recursive Green's function algorithm.
+"""
+import logging
+import numpy as np
+import nanonet.negf as negf
+import matplotlib.pyplot as plt
+from nanonet.tb import Hamiltonian
+from nanonet.tb import Orbitals
+from nanonet.tb.sorting_algorithms import sort_projection
+
+
+def main():
+    # 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"
+    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, 6, 22, 7, 23, 14, 30, 15, 31]
+
+    # define leads indices
+    hamiltonian = Hamiltonian(xyz=path,
+                              nn_distance=2.4,
+                              sort_func=sort_projection,
+                              left_lead=left_lead,
+                              right_lead=right_lead).initialize()
+
+    # set periodic boundary conditions
+    a_si = 5.50
+    primitive_cell = [[0, 0, a_si]]
+    hamiltonian.set_periodic_bc(primitive_cell)
+
+    # get Hamiltonian matrices
+    hl_bd, h0_bd, hr_bd, subblocks = hamiltonian.get_hamiltonians_block_tridiagonal(optimized=True)
+
+    # specify energy array
+    energy = np.linspace(2.1, 2.5, 50)
+
+    # specify dephasing constant
+    damp = 0.001j
+
+    # 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("{0} Energy: {1:.4f} eV".format(j, E))
+
+        # compute self-energies describing boundary conditions at the leads contacts
+        R, L = negf.surface_greens_function(E, hl_bd, h0_bd, hr_bd, iterate=True, damp=damp)
+
+        # compute Green's functions using the recursive Green's function algorithm
+        g_trans, grd, grl, gru, gr_left = negf.recursive_gf(E,
+                                                            hl_bd,
+                                                            h0_bd,
+                                                            hr_bd,
+                                                            left_se=L,
+                                                            right_se=R,
+                                                            damp=damp)
+        # number of subblocks
+        num_blocks = len(grd)
+
+        # compute DOS
+        for jj in range(num_blocks):
+            dos[j] = dos[j] - np.trace(np.imag(grd[jj])) / num_blocks
+
+        # coupling matrices
+        gamma_l = 1j * (L - L.conj().T)
+        gamma_r = 1j * (R - R.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
+    fig, ax = plt.subplots(2, 1)
+    ax[0].plot(energy, dos, 'k')
+    ax[0].set_ylabel(r'DOS (a.u)')
+    ax[0].set_xlabel(r'Energy (eV)')
+
+    ax[1].plot(energy, tr, 'k')
+    ax[1].set_ylabel(r'Transmission (a.u.)')
+    ax[1].set_xlabel(r'Energy (eV)')
+    fig.tight_layout()
+    plt.show()
+
+
+if __name__ == '__main__':
+
+    main()

nanonet/negf/greens_functions.py

@@ -3,19 +3,22 @@
 """
 from __future__ import print_function, division
 import numpy as np
+import scipy.sparse as sp
 import scipy.linalg as linalg
 
 
-def surface_greens_function_poles(h_list):
+def surface_greens_function_poles(hl, h0, hr):
     """Computes eigenvalues and eigenvectors for the complex band structure problem.
     The eigenvalues correspond to the wave vectors as `exp(ik)`.
 
     Parameters
     ----------
-    h_list : list
-        List of the Hamiltonian blocks - blocks describes coupling
-        with left-side lead, Hamiltonian of the device region and
-        coupling with right-side lead
+    hl : numpy.ndarray (dtype=numpy.float)
+        Left-side coupling Hamiltonian
+    h0 : numpy.ndarray (dtype=numpy.float)
+        Hamiltonian of the device region
+    hr : numpy.ndarray (dtype=numpy.float)
+        Right-side coupling Hamiltonian
 
     Returns
     -------
@@ -26,27 +29,48 @@ def surface_greens_function_poles(h_list):
     """
 
     # linearize polynomial eigenvalue problem
-    pr_order = len(h_list) - 1
-    matix_size = h_list[0].shape[0]
-    full_matrix_size = pr_order * matix_size
-    identity = np.identity(matix_size)
 
+    shapes_d = []
+    shapes_ud = []
+    shapes_ld = []
+
+    if isinstance(h0, list):
+        num_blocks = len(h0)
+        for j in range(len(h0)):
+            shapes_d.append(h0[j].shape[0])
+            shapes_ld.append(hl[j].shape)
+            shapes_ud.append(hr[j].shape)
+        matix_size = np.sum(shapes_d)
+    else:
+        num_blocks = 1
+        shapes_d.append(h0.shape[0])
+        shapes_ld.append(hl.shape)
+        shapes_ud.append(hr.shape)
+        matix_size = h0.shape[0]
+        hl = [hl]
+        h0 = [h0]
+        hr = [hr]
+
+    full_matrix_size = 2 * matix_size
+    identity = np.identity(matix_size)
     main_matrix = np.zeros((full_matrix_size, full_matrix_size), dtype=np.complex)
     overlap_matrix = np.zeros((full_matrix_size, full_matrix_size), dtype=np.complex)
+    main_matrix[0:matix_size, matix_size:2 * matix_size] = identity
+    overlap_matrix[0:matix_size, 0:matix_size] = identity
+    main_matrix[matix_size:matix_size+hl[-1].shape[0], matix_size-hl[-1].shape[1]:matix_size] = -hl[-1]
+    overlap_matrix[2*matix_size-hr[-1].shape[0]:2*matix_size, matix_size:matix_size+hr[-1].shape[1]] = hr[-1]
+
+    offfset = matix_size
+
+    for j in range(num_blocks):
 
-    for j in range(pr_order):
+        main_matrix[offfset:offfset+h0[j].shape[0], offfset:offfset+h0[j].shape[0]] = -h0[j]
 
-        main_matrix[(pr_order - 1) * matix_size:pr_order * matix_size,
-                    j * matix_size:(j + 1) * matix_size] = -h_list[j]
+        if j > 0:
+            main_matrix[offfset:offfset + hl[j - 1].shape[0], offfset - hl[j - 1].shape[1]:offfset] = hl[j - 1]
+            main_matrix[offfset - hl[j - 1].shape[1]:offfset, offfset:offfset + hl[j - 1].shape[0]] = hr[j - 1]
 
-        if j == pr_order - 1:
-            overlap_matrix[j * matix_size:(j + 1) * matix_size,
-                           j * matix_size:(j + 1) * matix_size] = h_list[pr_order]
-        else:
-            overlap_matrix[j * matix_size:(j + 1) * matix_size,
-                           j * matix_size:(j + 1) * matix_size] = identity
-            main_matrix[j * matix_size:(j + 1) * matix_size,
-                        (j + 1) * matix_size:(j + 2) * matix_size] = identity
+        offfset += h0[j].shape[0]
 
     alpha, betha, _, eigenvects, _, _ = linalg.lapack.cggev(main_matrix, overlap_matrix)
 
@@ -116,8 +140,7 @@ def iterate_gf(E, h_0, h_l, h_r, se, num_iter):
     """
 
     for _ in range(num_iter):
-        se = h_r.dot(np.linalg.pinv(E * np.identity(h_0.shape[0]) - h_0 - se)).dot(h_l)
-
+        se = h_r.dot(linalg.pinv(E * np.identity(h_0.shape[0]) - h_0 - se)).dot(h_l)
     return se
 
 
@@ -150,19 +173,87 @@ def surface_greens_function(E, h_l, h_0, h_r, iterate=False, damp=0.0001j):
         Right-lead self-energy matrix
     """
 
-    h_list = [h_l, h_0 - (E+damp) * np.identity(h_0.shape[0]), h_r]
-    vals, vects = surface_greens_function_poles(h_list)
+    # linearize polynomial eigenvalue problem
+
+    shapes_d = []
+    shapes_ud = []
+    shapes_ld = []
+
+    if isinstance(h_0, list):
+        num_blocks = len(h_0)
+        for j in range(len(h_0)):
+            shapes_d.append(h_0[j].shape[0])
+            shapes_ld.append(h_l[j].shape)
+            shapes_ud.append(h_r[j].shape)
+        matix_size = np.sum(shapes_d)
+    else:
+        num_blocks = 1
+        shapes_d.append(h_0.shape[0])
+        shapes_ld.append(h_l.shape)
+        shapes_ud.append(h_r.shape)
+        matix_size = h_0.shape[0]
+        h_l = [h_l]
+        h_0 = [h_0]
+        h_r = [h_r]
+
+    full_matrix_size = 2 * matix_size
+    identity = np.identity(matix_size)
+    main_matrix = np.zeros((full_matrix_size, full_matrix_size), dtype=np.complex)
+    overlap_matrix = np.zeros((full_matrix_size, full_matrix_size), dtype=np.complex)
+    main_matrix[0:matix_size, matix_size:2 * matix_size] = identity
+    overlap_matrix[0:matix_size, 0:matix_size] = identity
+    main_matrix[matix_size:matix_size+h_l[-1].shape[0], matix_size-h_l[-1].shape[1]:matix_size] = -h_l[-1]
+    overlap_matrix[2*matix_size-h_r[-1].shape[0]:2*matix_size, matix_size:matix_size+h_r[-1].shape[1]] = h_r[-1]
+
+    offfset = matix_size
+
+    for j in range(num_blocks):
+
+        main_matrix[offfset:offfset+h_0[j].shape[0], offfset:offfset+h_0[j].shape[0]] = -h_0[j]
+
+        if j > 0:
+            main_matrix[offfset:offfset + h_l[j - 1].shape[0], offfset - h_l[j - 1].shape[1]:offfset] = -h_l[j - 1]
+            main_matrix[offfset - h_l[j - 1].shape[1]:offfset, offfset:offfset + h_l[j - 1].shape[0]] = -h_r[j - 1]
+
+        offfset += h_0[j].shape[0]
+
+    main_matrix[matix_size:2*matix_size, matix_size:2*matix_size] = \
+        (E + damp) * np.identity(matix_size) + main_matrix[matix_size:2*matix_size, matix_size:2*matix_size]
+
+    alpha, betha, _, vects, _, _ = linalg.lapack.cggev(main_matrix, overlap_matrix)
+
+    betha[betha == 0] = 1e-19
+    vals = alpha / betha
+
+    # sort absolute values
+    ind = np.argsort(np.abs(vals))
+    vals = vals[ind]
+    vects = vects[:, ind]
+    vects = vects[matix_size:, :]
+
+    # vals, vects = surface_greens_function_poles(E, h_l, h_0, h_r, damp)
     num_eigs = len(vals)
 
-    u_right = vects[:, :num_eigs//2]
-    u_left = vects[:, num_eigs-1:num_eigs//2-1:-1]
-    lambda_right = np.diag(vals[:num_eigs//2])
-    lambda_left = np.diag(vals[num_eigs-1:num_eigs//2-1:-1])
+    u_right = vects[:, :num_eigs // 2]
+    u_left = vects[:, num_eigs - 1:num_eigs // 2 - 1:-1]
+    lambda_right = np.diag(vals[:num_eigs // 2])
+    lambda_left = np.diag(vals[num_eigs - 1:num_eigs // 2 - 1:-1])
 
-    sgf_l = h_r.dot(u_right).dot(lambda_right).dot(np.linalg.pinv(u_right))
-    sgf_r = h_l.dot(u_left).dot(np.linalg.inv(lambda_left)).dot(np.linalg.pinv(u_left))
+    h0 = -main_matrix[matix_size:2*matix_size, matix_size:2*matix_size] + (E + damp) * np.identity(matix_size)
+    hl = -main_matrix[matix_size:2*matix_size, 0:matix_size]
+    hr = overlap_matrix[matix_size:2*matix_size, matix_size:2*matix_size]
+
+    sgf_l = hr.dot(u_right).dot(lambda_right).dot(np.linalg.pinv(u_right))
+    sgf_r = hl.dot(u_left).dot(np.linalg.inv(lambda_left)).dot(np.linalg.pinv(u_left))
 
     if iterate:
-        return iterate_gf(E, h_0, h_l, h_r, sgf_l, 2), iterate_gf(E, h_0, h_r, h_l, sgf_r, 2)
+        sgf_l = iterate_gf(E, h0, hl, hr, sgf_l, 2)
+        sgf_r = iterate_gf(E, h0, hr, hl, sgf_r, 2)
+
+    s01, s02 = h_0[0].shape
+    s11, s12 = h_0[-1].shape
+
+    sgf_l = sgf_l[-s11:, -s12:]
+    sgf_r = sgf_r[:s01, :s02]
 
     return sgf_l, sgf_r

nanonet/negf/recursive_greens_functions.py

@@ -27,7 +27,7 @@ def mat_left_div(mat_a, mat_b):
     return ans
 
 
-def recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=0, s_out=0, damp=0.000001j):
+def _recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=0, s_out=0, damp=0.000001j):
     """The recursive Green's function algorithm is taken from
     M. P. Anantram, M. S. Lundstrom and D. E. Nikonov, Proceedings of the IEEE, 96, 1511 - 1550 (2008)
     DOI: 10.1109/JPROC.2008.927355
@@ -204,3 +204,79 @@ def recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=0, s_out=0, da
                grd, grl, gru, gr_left, \
                gnd, gnl, gnu, gin_left, \
                gpd, gpl, gpu, gip_left
+
+
+def recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, left_se=None, right_se=None, s_in=0, s_out=0, damp=0.000001j):
+    """The recursive Green's function algorithm is taken from
+    M. P. Anantram, M. S. Lundstrom and D. E. Nikonov, Proceedings of the IEEE, 96, 1511 - 1550 (2008)
+    DOI: 10.1109/JPROC.2008.927355
+
+    In order to get the electron correlation function output, the parameters s_in has to be set.
+    For the hole correlation function, the parameter s_out has to be set.
+
+    Parameters
+    ----------
+    energy : numpy.ndarray (dtype=numpy.float)
+        Energy array
+    mat_d_list : list of numpy.ndarray (dtype=numpy.float)
+        List of diagonal blocks
+    mat_u_list : list of numpy.ndarray (dtype=numpy.float)
+        List of upper-diagonal blocks
+    mat_l_list : list of numpy.ndarray (dtype=numpy.float)
+        List of lower-diagonal blocks
+    s_in :
+         (Default value = 0)
+    s_out :
+         (Default value = 0)
+    damp :
+         (Default value = 0.000001j)
+
+    Returns
+    -------
+    g_trans : numpy.ndarray (dtype=numpy.complex)
+        Blocks of the retarded Green's function responsible for transmission
+    grd : numpy.ndarray (dtype=numpy.complex)
+        Diagonal blocks of the retarded Green's function
+    grl : numpy.ndarray (dtype=numpy.complex)
+        Lower diagonal blocks of the retarded Green's function
+    gru : numpy.ndarray (dtype=numpy.complex)
+        Upper diagonal blocks of the retarded Green's function
+    gr_left : numpy.ndarray (dtype=numpy.complex)
+        Left-conencted blocks of the retarded Green's function
+    gnd : numpy.ndarray (dtype=numpy.complex)
+        Diagonal blocks of the retarded Green's function
+    gnl : numpy.ndarray (dtype=numpy.complex)
+        Lower diagonal blocks of the retarded Green's function
+    gnu : numpy.ndarray (dtype=numpy.complex)
+        Upper diagonal blocks of the retarded Green's function
+    gin_left : numpy.ndarray (dtype=numpy.complex)
+        Left-conencted blocks of the retarded Green's function
+    gpd : numpy.ndarray (dtype=numpy.complex)
+        Diagonal blocks of the retarded Green's function
+    gpl : numpy.ndarray (dtype=numpy.complex)
+        Lower diagonal blocks of the retarded Green's function
+    gpu : numpy.ndarray (dtype=numpy.complex)
+        Upper diagonal blocks of the retarded Green's function
+    gip_left : numpy.ndarray (dtype=numpy.complex)
+        Left-conencted blocks of the retarded Green's function
+    """
+
+    if isinstance(left_se, np.ndarray):
+        s01, s02 = mat_d_list[0].shape
+        left_se = left_se[:s01, :s02]
+        mat_d_list[0] = mat_d_list[0] + left_se
+
+    if isinstance(right_se, np.ndarray):
+        s11, s12 = mat_d_list[-1].shape
+        right_se = right_se[-s11:, -s12:]
+        mat_d_list[-1] = mat_d_list[-1] + right_se
+
+    ans = _recursive_gf(energy, mat_l_list, mat_d_list, mat_u_list, s_in=s_in, s_out=s_out, damp=damp)
+
+    if isinstance(left_se, np.ndarray):
+        mat_d_list[0] = mat_d_list[0] - left_se
+
+    if isinstance(right_se, np.ndarray):
+        mat_d_list[-1] = mat_d_list[-1] - right_se
+
+    return ans