Examine LDOS and bad band structure of spatial kwant systems

code/modules/analysis_funcs.py

@@ -68,7 +68,7 @@ def plot_bulk_gap(system_params, eigval_number=0, xlims=None, ylims=None):
     print(np.min(np.abs(eigvals)))
     
     
-def ldos_from_kpm(fsyst, params, energy=0, fig_size=(8, 6), vmax=1):
+def ldos_from_kpm(fsyst, params, energy=0, fig_size=(8, 6), vmax=1, norbs=6):
     """
     Plots the local density of states of a given system using KPM methods.
     

code/modules/potentials.py

@@ -287,7 +287,7 @@ def potential_system(L,W,lattice_spacing,params):
     :param list of dicts potential: patch or chain to be created 
     '''
 
-    lat = kwant.lattice.square(lattice_spacing) #builds the lattice
+    lat = kwant.lattice.square(lattice_spacing,norbs=1) #builds the lattice
     syst = kwant.Builder() #initializes the Hamiltonian
     
 
@@ -297,25 +297,24 @@ def potential_system(L,W,lattice_spacing,params):
 
     return syst
 
-def complex_potential_system(L,W,lattice_spacing,params):
+def infinite_potential_system(L,W,lattice_spacing,params):
     '''
-    Create a kwant system that includes two helicity sublattices of the 2D p orbitals, and distinguishes different hoppings between the different elements of the trimers.
+    :param float L: length of substrate
+    :param float W: width of substrate
+    :param float lattice_spacing: lattice spacing of the substrate
+    :param list of dicts potential: patch or chain to be created 
     '''
     
-    lat = kwant.lattice.Polyatomic(prim_vecs = [[1,0],[0,1]], basis = [[0,0],[0,0], [0,0],[0,0], [0,0],[0,0]], norbs = 1)
-    a_lat_pos, a_lat_neg, b_lat_pos, b_lat_neg, c_lat_pos, c_lat_neg = lat.sublattices
-    syst = kwant.Builder() #initializes the Hamiltonian
-    
-    complex_phase = np.exp(1j*2*params['phi'])
-    
-    #give all the sublattices the correct shape
-    for sublattice in lat.sublattices:
-        syst[(sublattice(i, j) for i in range(L) for j in np.linspace(0,-W,W,endpoint=False))] = onsite
-        syst[sublattice.neighbors()] = -t #off diagonal elements to the hopping
+    syst = kwant.Builder(symmetry=kwant.lattice.TranslationalSymmetry([1,-1]))
+
+    lat = kwant.lattice.square(lattice_spacing) #builds the lattice
     
+
+    syst[(lat(i, j) for i in range(L) for j in np.linspace(0,-W,W,endpoint=False))] = onsite #equates the diagonal elements of the Hamiltonian to the potential
+    syst[lat.neighbors()] = -t #off diagonal elements to the hopping
     syst = syst.finalized() #hamiltonian is built
-    return syst
 
+    return syst
 
 def sorted_evs(syst, potential, k=15):
     '''
@@ -368,10 +367,15 @@ def plotter(X,Y,sorted_evecs,sorted_evals,L,W, aspect = None,rounding=2):
             if aspect:
                 axs.set_aspect(aspect)
     
-def view_potential_shape(potential_type, n=2,Nc=1,a=6.75,Vo=0.306,closed_chain=False,L=100,W=100, inversion_symmetric=True,trivial=False):
+def view_potential_shape(potential_type, n=2,Nc=1,a=6.75,Vo=0.306,closed_chain=False,L=100,W=100, inversion_symmetric=True,trivial=False,view=True,semi_infinite=False):
+    
+    if semi_infinite:
+        syst = kwant.Builder(symmetry=kwant.lattice.TranslationalSymmetry([1,-1]))
+    else:
+        syst = kwant.Builder()
+        
+    lat = kwant.lattice.square(1, norbs=1)
     
-    lat = kwant.lattice.square(1)
-    syst = kwant.Builder()
 
     def shape_func(pos):
         x,y = pos
@@ -390,8 +394,11 @@ def view_potential_shape(potential_type, n=2,Nc=1,a=6.75,Vo=0.306,closed_chain=F
             else:
                 return False
 
-    syst[lat.shape(shape_func, (0,0) )] = 1
-    syst[lat.neighbors(n=1)] = 1
-    kwant.plot(syst)
+    syst[lat.shape(shape_func, (0,0) )] = 0
+    syst[lat.neighbors(n=1)] = -t
+    if view == True:
+        kwant.plot(syst)
+    else:
+        return syst
 
 

code/modules/potentials_rings.py

@@ -122,7 +122,10 @@ def chain_n(site, n, Nc, a, Vo, start_pos, closed_chain, inversion_symmetric = T
     :param bool inversion_symmetric: if False, then connectors 1-2 have a different size than connectors 1-3 in trimers_n.
     '''
     #connecting patch size
-    n_ = 1#n-1-n%2
+    long_dim = (n-1-n%2)/2
+    short_dim = long_dim/3*2
+    
+    n_ = long_dim
     
     for Uc in range(Nc):
         if Uc != 0:
@@ -138,10 +141,10 @@ def chain_n(site, n, Nc, a, Vo, start_pos, closed_chain, inversion_symmetric = T
         extra_ring = ring_n(site,a, Vo, n, [start_pos[0]+Nc*a*(n+n_), start_pos[1]+Nc*a*(n+n_)])
         if extra_ring < Vo:
             return extra_ring
-        extra_connector_1_2 = patch_n(site,a, Vo, n_, [start_pos[0]+(Nc*(n+n_)+1)*a, start_pos[1]+(Nc*(n+n_)-n_)*a])
+        extra_connector_1_2 = patch_n(site,a, Vo, short_dim, [start_pos[0]+(Nc*(n+n_)+1)*a, start_pos[1]+(Nc*(n+n_)-n_)*a],ydim=long_dim*a)
         if extra_connector_1_2 < Vo:
             return extra_connector_1_2
-        extra_connector_1_3 = patch_n(site,a, Vo, n_, [start_pos[0]+(Nc*(n+n_)-n_)*a, start_pos[1]+(Nc*(n+n_)+1)*a])
+        extra_connector_1_3 = patch_n(site,a, Vo, long_dim, [start_pos[0]+(Nc*(n+n_)-n_)*a, start_pos[1]+(Nc*(n+n_)+1)*a],ydim=short_dim*a)
         if extra_connector_1_3 < Vo:
             return extra_connector_1_3
             
@@ -168,8 +171,8 @@ def trimer_n(site,a,Vo,n,start_pos, in_chain_bulk=False, inversion_symmetric = T
     assert n >= 2
     
     #connecting patch size
-    long_dim = 1#n-1-n%2
-    short_dim = 0.25
+    long_dim = (n-1-n%2)/2
+    short_dim = long_dim/3*2
     
     n_ = long_dim
     
@@ -193,33 +196,33 @@ def trimer_n(site,a,Vo,n,start_pos, in_chain_bulk=False, inversion_symmetric = T
         return 0
     
     if n%2 == 0:
-        patch_1_2 = patch_n(site,a, Vo, long_dim, [start_pos[0]+n*a, start_pos[1]],ydim = short_dim*a)
+        patch_1_2 = patch_n(site,a, Vo, long_dim, [start_pos[0]+n*a, start_pos[1]+(n/2-short_dim/2)*a],ydim = short_dim*a)
         if patch_1_2 == 0:
             return 0
-        patch_1_3 = patch_n(site,a, Vo, short_dim, [start_pos[0], start_pos[1]+n*a],ydim=long_dim*a)
+        patch_1_3 = patch_n(site,a, Vo, short_dim, [start_pos[0]+(n/2-short_dim/2)*a, start_pos[1]+n*a],ydim=long_dim*a)
         if patch_1_3 == 0:
             return 0
     else:
-        patch_1_2 = patch_n(site,a, Vo, long_dim, [start_pos[0]+n*a, start_pos[1]+1*a],ydim=short_dim*a)
+        patch_1_2 = patch_n(site,a, Vo, long_dim, [start_pos[0]+n*a, start_pos[1]+(n/2-short_dim/2)*a],ydim=short_dim*a)
         if patch_1_2 == 0:
             return 0
-        patch_1_3 = patch_n(site,a, Vo, short_dim, [start_pos[0]+1*a, start_pos[1]+n*a],ydim=long_dim*a)
+        patch_1_3 = patch_n(site,a, Vo, short_dim, [start_pos[0]+(n/2-short_dim/2)*a, start_pos[1]+n*a],ydim=long_dim*a)
         if patch_1_3 == 0:
             return 0
     
     if in_chain_bulk:
         if n%2 == 0:
-            extra_connector_1_2 = patch_n(site,a, Vo, short_dim, [start_pos[0]+1*a, start_pos[1]-n_*a],ydim=long_dim*a)
+            extra_connector_1_2 = patch_n(site,a, Vo, short_dim, [start_pos[0]+(n/2-short_dim/2)*a, start_pos[1]-n_*a],ydim=long_dim*a)
             if extra_connector_1_2 == 0:
                 return 0
-            extra_connector_1_3 = patch_n(site,a, Vo, long_dim, [start_pos[0]-n_*a, start_pos[1]+1*a],ydim=short_dim*a)
+            extra_connector_1_3 = patch_n(site,a, Vo, long_dim, [start_pos[0]-n_*a, start_pos[1]+(n/2-short_dim/2)*a],ydim=short_dim*a)
             if extra_connector_1_3 == 0:
                 return 0
         else:
-            extra_connector_1_2 = patch_n(site,a, Vo, short_dim, [start_pos[0]+1*a, start_pos[1]-n_*a],ydim=long_dim*a)
+            extra_connector_1_2 = patch_n(site,a, Vo, short_dim, [start_pos[0]+(n/2-short_dim/2)*a, start_pos[1]-n_*a],ydim=long_dim*a)
             if extra_connector_1_2 == 0:
                 return 0
-            extra_connector_1_3 = patch_n(site,a, Vo, long_dim, [start_pos[0]-n_*a, start_pos[1]+1*a],ydim=short_dim*a)
+            extra_connector_1_3 = patch_n(site,a, Vo, long_dim, [start_pos[0]-n_*a, start_pos[1]+(n/2-short_dim/2)*a],ydim=short_dim*a)
             if extra_connector_1_3 == 0:
                 return 0
     return Vo
@@ -272,7 +275,7 @@ def ring_n(site, a, Vo, n, start_pos):
     r = np.sqrt(x**2+y**2)
     
     
-    if r <= n*a/2 and r>= (n-n/2)*a/2:
+    if r <= n*a/2 and r>= (n-n/4)*a/2:
         return 0
     else:
         return Vo
@@ -313,6 +316,25 @@ def potential_system(L,W,lattice_spacing,params):
 
     return syst
 
+def infinite_potential_system(L,W,lattice_spacing,params):
+    '''
+    :param float L: length of substrate
+    :param float W: width of substrate
+    :param float lattice_spacing: lattice spacing of the substrate
+    :param list of dicts potential: patch or chain to be created 
+    '''
+    
+    syst = kwant.Builder(symmetry=kwant.lattice.TranslationalSymmetry([1,-1]))
+
+    lat = kwant.lattice.square(lattice_spacing) #builds the lattice
+    
+
+    syst[(lat(i, j) for i in range(L) for j in np.linspace(0,-W,W,endpoint=False))] = onsite #equates the diagonal elements of the Hamiltonian to the potential
+    syst[lat.neighbors()] = -t #off diagonal elements to the hopping
+    syst = syst.finalized() #hamiltonian is built
+
+    return syst
+
 
 
 
@@ -391,8 +413,44 @@ def view_potential_shape(potential_type, n=2,Nc=1,a=6.75,Vo=0.306,closed_chain=F
             else:
                 return False
 
-    syst[lat.shape(shape_func, (int((n-n/4)*a),-int((n-n/4)*a)) )] = 1
+    syst[lat.shape(shape_func, (int((n-n/8)*a),-int((n-n/8)*a)) )] = 1
     syst[lat.neighbors(n=1)] = 1
     kwant.plot(syst)
 
 
+def view_potential_shape(potential_type, n=2,Nc=1,a=6.75,Vo=0.306,closed_chain=False,L=100,W=100, inversion_symmetric=True,trivial=False,view=True,semi_infinite=False):
+    
+    if semi_infinite:
+        syst = kwant.Builder(symmetry=kwant.lattice.TranslationalSymmetry([1,-1]))
+    else:
+        syst = kwant.Builder()
+        
+    lat = kwant.lattice.square(1, norbs=1)
+    
+
+    def shape_func(pos):
+        x,y = pos
+        if x<0 or x> L or y>0 or y<-W:
+            return False
+        else:
+            if potential_type == 'trimer':
+                res = trimer_n(pos,a=a,Vo=Vo,n=n,start_pos=[0,0], in_chain_bulk=False,inversion_symmetric=inversion_symmetric)
+                
+            elif potential_type == 'chain':
+                res = chain_n(pos, n=n, Nc=Nc, a=a, Vo=Vo, start_pos=[0,0], closed_chain=closed_chain,inversion_symmetric=inversion_symmetric)
+            elif potential_type == 'ssh':
+                res = ssh_n(pos, n=n, a=a, Vo=Vo, start_pos=[0,0], trivial=trivial)
+            elif potential_type == 'ring':
+                res = ring_n(pos, a=a, Vo=Vo, n=n, start_pos=[0,0])
+            if res == 0:
+                    return True
+            else:
+                return False
+
+    syst[lat.shape(shape_func, (int((n-n/8)*a),-int((n-n/8)*a)) )] = 0
+    syst[lat.neighbors(n=1)] = -t
+    if view == True:
+        kwant.plot(syst)
+    else:
+        return syst
+

code/potentials_chain.py

+# +
 import kwant
 import proplot as plot
 import scipy.sparse.linalg as sla
 import numpy as np
 import modules.potentials as potentials
 import modules.potentials_rings as potentials_rings
+
+from modules.analysis_funcs import ldos_from_kpm
 from matplotlib import pyplot as plt
 
-potentials_rings.view_potential_shape('chain',Nc=5,closed_chain=True,n=3,L=100,W=100, inversion_symmetric=False,trivial=True)
+# +
+# potentials.view_potential_shape('chain',Nc=3,closed_chain=False,n=3,L=200,W=200, inversion_symmetric=False,trivial=True)
+
+# +
+hbar = 1
+m = 1
+a = 6.75 #0.36 nm (lattice constant) in a.u.
+
+
+lattice_spacing = 1
+size = 3
+Nc = 3
+
+
+size_ = size-1-size%2
+L = int( (size+size_)*a*Nc + size*a )
+W = L#sets the size of the discretized lattice (number points)
+chain = dict(
+    L=L,
+    W=W,
+    size=3,
+    a = a,
+    start_pos=[0,0],
+    Nc = 0,
+    Vo = 0.5,
+    closed_chain=False,
+    trivial=False,
+    inversion_symmetric=True,
+    potential_type = 'chain_n'
+)
+# -
+
+inf_chain = potentials.view_potential_shape(L=L,W=W,n=3,Nc=3,a=a,Vo=0.5,view=False,potential_type='chain',semi_infinite=True)
+
+p = kwant.plotter.bands(inf_chain.finalized(), momenta = np.linspace(-np.pi,np.pi,1000),params=chain)
+
+fin_chain = potentials_rings.view_potential_shape(L=150,W=150,n=3,Nc=5,a=a,Vo=0.5,view=False,potential_type='chain',semi_infinite=False)
+
+ldos_from_kpm(fin_chain.finalized(), chain, vmax = 3, energy = 1.75)
 
 # # Patches
 
@@ -24,7 +65,7 @@ well_pos = 2*a #sets the spacing between the wells
 chain = dict(
     L=L,
     W=W,
-    size=5,
+    size=3,
     start_pos = [0,0],
     Nc = 0,
     closed_chain = False,
@@ -53,17 +94,21 @@ a = 6.75 #0.36 nm (lattice constant) in a.u.
 
 
 lattice_spacing = 1
-L, W = 100, 100 #sets the size of the discretized lattice (number points)
-well_pos = 2*a #sets the spacing between the wells
+size = 3
+Nc = 1
 
+
+size_ = size-1-size%2
+L = int( (size+size_)*a*Nc + size*a )
+W = L#sets the size of the discretized lattice (number points)
 chain = dict(
     L=L,
     W=W,
-    size=5,
+    size=3,
     a = a,
     start_pos=[0,0],
     Nc = 0,
-    Vo = 0.0506164,
+    Vo = 0.506164,
     closed_chain=False,
     trivial=False,
     inversion_symmetric=True,
@@ -71,7 +116,7 @@ chain = dict(
 )
 
 # +
-syst = potentials_rings.potential_system(L,W,lattice_spacing,chain)
+syst = potentials.potential_system(L,W,lattice_spacing,chain)
 sorted_evals, sorted_evecs = potentials.sorted_evs(syst, chain,k=30)
 
 x, y = np.linspace(0,L,L,endpoint=False), np.linspace(0, -W, W, endpoint=False)
@@ -88,7 +133,7 @@ a = 6.75 #0.36 nm (lattice constant) in a.u.
 lattice_spacing = 1
 
 size = 3
-Nc = 3
+Nc = 7
 
 
 size_ = size-1-size%2
@@ -104,7 +149,7 @@ chain = dict(
     size=size,
     start_pos=[0,0],
     Nc=Nc,
-    Vo = 0.0506164,
+    Vo = 0.506164,
     closed_chain=False,
     inversion_symmetric=True,
     potential_type = 'chain_n',
@@ -112,7 +157,7 @@ chain = dict(
 )
 
 # +
-syst = potentials_rings.potential_system(L,W,lattice_spacing,chain)
+syst = potentials.potential_system(L,W,lattice_spacing,chain)
 sorted_evals, sorted_evecs = potentials.sorted_evs(syst, chain,k=100)
 
 x, y = np.linspace(0,L,L,endpoint=False), np.linspace(0, -W, W, endpoint=False)