Merge pull request #77 from freude/develop

corrected LDOS in bi_nanoribbon_transport.py

examples/bi_nanoribbon_projected.py

@@ -54,7 +54,7 @@ def main():
     #          np.diag([1, 1, 1, 1] * h.num_of_nodes, -4))
     # proj_y = np.matrix(np.diag([-1j, -1j, -1j, -1j] * h.num_of_nodes, 4) + \
     #          np.diag([1j, 1j, 1j, 1j] * h.num_of_nodes, 4))
-    proj_z = np.matrix(np.diag([1, 1, 1, 1, -1, -1, -1, -1] * h.num_of_nodes))
+    proj_z = np.matrix(np.diag([0, 1, 1, 1, 0, -1, -1, -1] * h.num_of_nodes))
     proj_s = np.matrix(np.diag([1, 0, 0, 0, 1, 0, 0, 0] * h.num_of_nodes))
 
     # --------------------------------------------------------------
@@ -95,7 +95,7 @@ def main():
 
         im1 = axs[0].scatter(k_points[:, 0], item, marker='o',
                              s=100*spin_structure1[:, j], facecolors='none', edgecolors='r')
-        im1 = axs[0].scatter(k_points[:, 0], item, marker='v',
+        im2 = axs[0].scatter(k_points[:, 0], item, marker='v',
                              s=100*spin_structure2[:, j], facecolors='none', edgecolors='r')
     axs[0].set_title('Magnitude of the z-axis spin projection', fontsize=10)
     axs[0].set_ylim((-1, 1))

examples/bi_nanoribbon_transport.py

@@ -22,7 +22,7 @@ def fd(energy, ef, temp):
         kb = 8.61733e-5  # Boltzmann constant in eV
         return 1.0 / (1.0 + np.exp((energy - ef) / (kb * temp)))
 
-    path_to_xyz_file = 'input_samples/bi_nanoribbon_090.xyz'
+    path_to_xyz_file = 'input_samples/bi_nanoribbon_014.xyz'
 
     bi = tb.Orbitals('Bi')
     bi.add_orbital("s", energy=-10.906, principal=0, orbital=0, magnetic= 0, spin=0)
@@ -77,7 +77,7 @@ def fd(energy, ef, temp):
 
         ind = np.argsort(np.array(h._coords)[:, 1])
         # gn_diag = np.concatenate((np.diag(gnd[0])[ind], np.diag(gnd[1])[ind], np.diag(gnd[2])[ind]))
-        gn_diag = np.diag(gnd[0])
+        gn_diag = np.diag(np.imag(grd[0]))
         gn_diag = np.reshape(gn_diag, (h.num_of_nodes, -1))
         dens[j, :] = 2 * np.sum(gn_diag, axis=1)
 
@@ -103,6 +103,13 @@ def fd(energy, ef, temp):
     ax.set_xlabel(r'Energy (eV)')
     plt.show()
 
+    plt.figure(figsize=[5, 10])
+    ax = plt.axes()
+    ax.contourf(np.arange(14), energy, dens, 200, cmap='terrain')
+    ax.set_ylabel(r'Energy (eV)')
+    ax.set_xlabel(r'Site number')
+    plt.show()
+
     data_to_write = np.c_[energy[:, None], tr]
     np.savetxt(path_to_dat_file, np.c_[data_to_write])
 

examples/bi_nanoribbon_transport_scf.py

@@ -3,6 +3,7 @@
 from negf.recursive_greens_functions import recursive_gf
 from examples import data_bi_nanoribbon
 
+
 def radial_dep(coords):
 
     norm_of_coords = np.linalg.norm(coords)
@@ -16,7 +17,7 @@ def radial_dep(coords):
         return 100
 
 
-def test_gf(energy):
+def make_tb_matrices():
 
     path_to_xyz_file = 'input_samples/bi_nanoribbon_014.xyz'
 
@@ -42,10 +43,15 @@ def test_gf(energy):
     h.set_periodic_bc([period])
     h_l, h_0, h_r = h.get_coupling_hamiltonians()
 
-    diags = np.zeros((energy.shape[0], h.num_of_nodes), dtype=np.complex)
+    return h_l, h_0, h_r, h.num_of_nodes
+
+
+def test_gf(energy, h_l, h_0, h_r, num_of_nodes):
+
+    diags = np.zeros((energy.shape[0], num_of_nodes), dtype=np.complex)
 
     for j, E in enumerate(energy):
-        L, R = tb.surface_greens_function(E, h_l, h_0, h_r, iterate=5)
+        L, R = tb.surface_greens_function(E, h_l, h_0, h_r, iterate=2)
 
         g_trans, grd, grl, gru, gr_left = recursive_gf(E,
                                                        [h_l],
@@ -54,7 +60,7 @@ def test_gf(energy):
 
         gn_diag = np.diag(grd[0])
         # gn_diag = np.reshape(gn_diag, (h._orbitals_dict['Bi'].num_of_orbitals, -1))
-        gn_diag = np.reshape(gn_diag, (h.num_of_nodes, -1))
+        gn_diag = np.reshape(gn_diag, (num_of_nodes, -1))
         diags[j, :] = 2 * np.sum(gn_diag, axis=1)
 
     return diags
@@ -63,19 +69,21 @@ def test_gf(energy):
 if __name__ == '__main__':
 
     import matplotlib.pyplot as plt
-    from tb.cfr import CFR
+    from negf.cfr import CFR
+
+    h_r, h_0, h_l, num_of_nodes = make_tb_matrices()
 
     tempr = 300
     fd = CFR(40)
     Ef = np.linspace(-1.0, 1.0, 40)
     ans = []
-    moment = 1j * CFR.val_inf * test_gf(np.array([1j * CFR.val_inf]))
+    moment = 1j * CFR.val_inf * test_gf(np.array([1j * CFR.val_inf]), h_l, h_0, h_r, num_of_nodes)
 
     for ef in Ef:
         print(ef)
         points = fd.genetate_integration_points(ef, tempr)
-        gf_vals = test_gf(points)
+        gf_vals = test_gf(points, h_l, h_0, h_r, num_of_nodes)
         ans.append(fd.integrate1(gf_vals, tempr, zero_moment=moment))
 
-    plt.plot(np.array(ans))
+    plt.plot(np.squeeze(np.array(ans)))
     plt.show()

negf/surf_greens_function.py

@@ -108,10 +108,9 @@ def surface_greens_function_poles(E, h_l, h_0, h_r):
 
     alpha, betha, _, eigenvects, _, _ = linalg.lapack.cggev(main_matrix, overlap_matrix)
 
-    eigenvals = np.zeros(alpha.shape, dtype=np.complex128)
+    eigenvals = np.zeros(alpha.shape, dtype=np.complex)
 
     for j, item in enumerate(zip(alpha, betha)):
-
         if np.abs(item[1]) != 0.0:
             eigenvals[j] = item[0] / item[1]
         else:

tb/greens_function.py

@@ -46,7 +46,7 @@ def surface_greens_function_poles(h_list):
 
     alpha, betha, _, eigenvects, _, _ = linalg.lapack.cggev(main_matrix, overlap_matrix)
 
-    eigenvals = np.zeros(alpha.shape, dtype=np.complex128)
+    eigenvals = np.zeros(alpha.shape, dtype=np.complex)
 
     for j, item in enumerate(zip(alpha, betha)):