Rework flow in potentials by adding a patch function. Add functionality to...

Rework flow in potentials by adding a patch function. Add functionality to close chain for potentials.

code/kwant_chain.py

@@ -26,32 +26,33 @@ system_params = dict(
     phi_d = 0.
 )
 
-# system_params = dict(
-#     mu_a = 0.,
-#     mu_b = 0.,
-#     mu_c = 0.,
-#     j2 = 1,
-#     phi_d = 0.
-# )
+system_params = dict(
+    mu_a = 0.,
+    mu_b = 0.,
+    mu_c = 0.,
+    j2 = 1,
+    phi_d = 0.
+)
 # -
 
-N_uc = 10
-closed_chain = False
-inf_chain = tilted_diamond_chain_system(N_uc,system_params, semi_infinite = True, closed_chain=closed_chain)
-fin_chain = tilted_diamond_chain_system(N_uc,system_params, semi_infinite = False, closed_chain=closed_chain)
-# inf_chain = s_tilted_diamond_chain_system(N_uc,system_params, semi_infinite = True, closed_chain=closed_chain)
-# fin_chain = s_tilted_diamond_chain_system(N_uc,system_params, semi_infinite = False, closed_chain=closed_chain)
+l=1
+N_uc = 3
+closed_chain = True
+# inf_chain = tilted_diamond_chain_system(l,N_uc,system_params, semi_infinite = True, closed_chain=closed_chain)
+# fin_chain = tilted_diamond_chain_system(l,N_uc,system_params, semi_infinite = False, closed_chain=closed_chain)
+inf_chain = s_tilted_diamond_chain_system(N_uc,system_params, semi_infinite = True, closed_chain=closed_chain)
+fin_chain = s_tilted_diamond_chain_system(N_uc,system_params, semi_infinite = False, closed_chain=closed_chain)
 o = kwant.plot(fin_chain)
 
 # ## LDOS
 
 # ### snapshots
 
-o = ldos_from_kpm(fin_chain.finalized(), system_params, vmax = 10, energy = 0 )
+o = ldos_from_kpm(fin_chain.finalized(), system_params, vmax = 20, energy = 0 )
 
-o = ldos_from_kpm(fin_chain.finalized(), system_params, vmax = 10, energy = 2*system_params['j2'] )
+o = ldos_from_kpm(fin_chain.finalized(), system_params, vmax = 30, energy = 2*system_params['j2'] )
 
-o = ldos_from_kpm(fin_chain.finalized(), system_params, vmax = 100, energy = 2*np.sqrt(2)*system_params['j2'] )
+o = ldos_from_kpm(fin_chain.finalized(), system_params, vmax = 230, energy = 2*np.sqrt(2)*system_params['j2'] )
 
 # ## Spectrum
 

code/modules/.ipynb_checkpoints/diamond_chain-checkpoint.py

@@ -110,7 +110,7 @@ def diamond_chain_system(N_c, system_params, semi_infinite = False, leads = Fals
     
     return syst
 
-def tilted_diamond_chain_system(N_c, system_params, semi_infinite = False, leads = False, closed_chain = False):
+def tilted_diamond_chain_system(l, N_c, system_params, semi_infinite = False, leads = False, closed_chain = False):
     '''
     Create a diamond chain of trimer unit cells. Each atom of the trimer has two orbital angular momentum states, + and -
           
@@ -128,6 +128,7 @@ def tilted_diamond_chain_system(N_c, system_params, semi_infinite = False, leads
                          
     The Φ represents an out-of-plane magnetic field. The phase is added along the == bond in each unit cell 
     
+    :param int l: the orbital angular momentum number. p bands have l=1, and d have l=2.
     :param int N_c: number of unit cells to include in the cell
     :param dict system_params: parameters
     :param bool semi_infinite: whether to make a semi-infinite chain or not
@@ -136,6 +137,8 @@ def tilted_diamond_chain_system(N_c, system_params, semi_infinite = False, leads
     
     :rtype kwant.system.FiniteSystem:
     '''
+    if l == 0:
+        raise ValueError('l cannot be 0: use s_tilted_diamond_chain_system instead.')
     
     # make lattices and sublattices
     lat = kwant.lattice.Polyatomic(prim_vecs = [[1,-1],[1,1]], basis = [[0,0],[0,0], [0,1],[0,1], [-1,0],[-1,0]], norbs = 1)
@@ -183,7 +186,7 @@ def tilted_diamond_chain_system(N_c, system_params, semi_infinite = False, leads
                 
             
         else:
-            added_phase = np.exp(1j*2*system_params['phi_d'])
+            added_phase = np.exp(1j*2*system_params['phi_d']*l)
             
         
         if i < N_c - 1:

code/modules/diamond_chain.py

@@ -110,7 +110,7 @@ def diamond_chain_system(N_c, system_params, semi_infinite = False, leads = Fals
     
     return syst
 
-def tilted_diamond_chain_system(N_c, system_params, semi_infinite = False, leads = False, closed_chain = False):
+def tilted_diamond_chain_system(l, N_c, system_params, semi_infinite = False, leads = False, closed_chain = False):
     '''
     Create a diamond chain of trimer unit cells. Each atom of the trimer has two orbital angular momentum states, + and -
           
@@ -128,6 +128,7 @@ def tilted_diamond_chain_system(N_c, system_params, semi_infinite = False, leads
                          
     The Φ represents an out-of-plane magnetic field. The phase is added along the == bond in each unit cell 
     
+    :param int l: the orbital angular momentum number. p bands have l=1, and d have l=2.
     :param int N_c: number of unit cells to include in the cell
     :param dict system_params: parameters
     :param bool semi_infinite: whether to make a semi-infinite chain or not
@@ -136,6 +137,8 @@ def tilted_diamond_chain_system(N_c, system_params, semi_infinite = False, leads
     
     :rtype kwant.system.FiniteSystem:
     '''
+    if l == 0:
+        raise ValueError('l cannot be 0: use s_tilted_diamond_chain_system instead.')
     
     # make lattices and sublattices
     lat = kwant.lattice.Polyatomic(prim_vecs = [[1,-1],[1,1]], basis = [[0,0],[0,0], [0,1],[0,1], [-1,0],[-1,0]], norbs = 1)
@@ -183,7 +186,7 @@ def tilted_diamond_chain_system(N_c, system_params, semi_infinite = False, leads
                 
             
         else:
-            added_phase = np.exp(1j*2*system_params['phi_d'])
+            added_phase = np.exp(1j*2*system_params['phi_d']*l)
             
         
         if i < N_c - 1:

code/modules/potentials.py

@@ -21,8 +21,7 @@ t = hbar**2/(2*m*lattice_spacing**2) #sets the hopping parameter (the off-diagon
 # $-{\frac  {\hbar ^{2}}{2m}}{\frac  {d^{2}\psi }{dx^{2}}}+V(x)\psi =E\psi \quad (1)$
 
 # +
-
-def chain_n(site, n, Nc, a, Vo, start_pos):
+def chain_n(site, n, Nc, a, Vo, start_pos, closed_chain):
     '''
     Create a chain of Nc unit cells composed of trimers of size n.
     
@@ -43,10 +42,20 @@ def chain_n(site, n, Nc, a, Vo, start_pos):
             if cell == 0:
                 return 0
             
+    if closed_chain:
+        extra_patch = patch_n(site,a, Vo, n, [start_pos[0]+Nc*a*(n+n_), start_pos[1]+Nc*a*(n+n_)])
+        if extra_patch == 0:
+            return 0
+        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])
+        if extra_connector_1_2 == 0:
+            return 0
+        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])
+        if extra_connector_1_3 == 0:
+            return 0
+            
     return Vo
-    
-    
-    
+
+
 def trimer_n(site,a,Vo,n,start_pos, in_chain_bulk=False):
     '''
     Create a trimer potential profile of nxn patches. If n is odd, connecting patches are (n-2)x(n-2) and centered,
@@ -69,57 +78,79 @@ def trimer_n(site,a,Vo,n,start_pos, in_chain_bulk=False):
     '''
     assert n >= 2
     
-    x,y = site.pos
-    x += -start_pos[0]
-    y += start_pos[1]
-    
     #connecting patch size
     n_ = n-1-n%2
     
-    #patch 1
-    if (0 <= x <= n*a and 0 >= y >= -n*a):
+    patch_1 = patch_n(site,a, Vo, n, start_pos)
+    if patch_1 == 0:
         return 0
-    #patch 2
-    elif ((n+n_)*a <= x <= (2*n+n_)*a and 0 >= y >= -n*a):
+    
+    patch_2 = patch_n(site,a, Vo, n, [start_pos[0]+(n+n_)*a, start_pos[1]])
+    if patch_2 == 0:
         return 0
-    #patch 3
-    elif (0 <= x <= n*a and -(n+n_)*a >= y >= -(2*n+n_)*a):
+    
+    patch_3 = patch_n(site,a, Vo, n, [start_pos[0], start_pos[1]+(n+n_)*a])
+    if patch_3 == 0:
         return 0
     
-    # connector between 1 and 2
-    if n%2==0:
-        if (n*a <= x <= (n+n_)*a and 0 >= y >= -n_*a):
+    if n%2 == 0:
+        patch_1_2 = patch_n(site,a, Vo, n_, [start_pos[0]+n*a, start_pos[1]])
+        if patch_1_2 == 0:
             return 0
-        # connector between 1 and 3
-        if (0 <= x <= n_*a and -n*a >= y >= -(n+n_)*a):
+        patch_1_3 = patch_n(site,a, Vo, n_, [start_pos[0], start_pos[1]+n*a])
+        if patch_1_3 == 0:
             return 0
     else:
-        if (n*a <= x <= (n+n_)*a and -1*a >= y >= (-1-n_)*a):
+        patch_1_2 = patch_n(site,a, Vo, n_, [start_pos[0]+n*a, start_pos[1]+1*a])
+        if patch_1_2 == 0:
             return 0
-        # connector between 1 and 3
-        if 1*a <= x <= (n_+1)*a and -n*a >= y >= -(n+n_)*a:
+        patch_1_3 = patch_n(site,a, Vo, n_, [start_pos[0]+1*a, start_pos[1]+n*a])
+        if patch_1_3 == 0:
             return 0
-        
+    
     if in_chain_bulk:
-        if n%2==0:
-            if (-n_*a <= x <= 0 and -1*a >= y >= -n*a):
+        if n%2 == 0:
+            extra_connector_1_2 = patch_n(site,a, Vo, n_, [start_pos[0]+1*a, start_pos[1]-n_*a])
+            if extra_connector_1_2 == 0:
                 return 0
-            if (1*a <= x <= n*a and n_*a >= y >= 0):
+            extra_connector_1_3 = patch_n(site,a, Vo, n_, [start_pos[0]-n_*a, start_pos[1]+1*a])
+            if extra_connector_1_3 == 0:
                 return 0
         else:
-            if (-n_*a <= x <= 0 and -1*a >= y >= (n-1)*a):
+            extra_connector_1_2 = patch_n(site,a, Vo, n_, [start_pos[0]+1*a, start_pos[1]-n_*a])
+            if extra_connector_1_2 == 0:
                 return 0
-            if 1*a <= x <= (n-1)*a and n_*a >= y >= 0:
+            extra_connector_1_3 = patch_n(site,a, Vo, n_, [start_pos[0]-n_*a, start_pos[1]+1*a])
+            if extra_connector_1_3 == 0:
                 return 0
-        
     return Vo
+
+
+def patch_n(site, a, Vo, n, start_pos):
+    try:
+        x, y = site.pos
+        x_ = x + -start_pos[0]
+    except:
+        x,y = site 
+        
+    x += -start_pos[0]
+    y += start_pos[1]
+    
+    
+    if 0 <= x <= n*a and 0 >= y >= -n*a:
+        return 0
+    else:
+        return Vo
+    
     
     
-def onsite(site, Vo, potential_type, a, size, start_pos, Nc):
-    if potential_type == 'trimer_n':
+def onsite(site, Vo, potential_type, a, size, start_pos, Nc, closed_chain):
+    if potential_type == 'patch_n':
+        return 4*t+patch_n(site, a, Vo, size, start_pos)
+    elif potential_type == 'trimer_n':
         return 4*t + trimer_n(site, a, Vo,size, start_pos)
-    if potential_type == 'chain_n':
-        return 4*t + chain_n(site, size, Nc, a, Vo, start_pos)
+    elif potential_type == 'chain_n':
+        return 4*t + chain_n(site, size, Nc, a, Vo, start_pos, closed_chain)
     
 
 
@@ -138,7 +169,7 @@ def potential_system(L,W,lattice_spacing,params):
     
 
     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[lat.neighbors()] = -1j*t #off diagonal elements to the hopping
     syst = syst.finalized() #hamiltonian is built
 
     return syst
@@ -177,5 +208,29 @@ def plotter(X,Y,sorted_evecs,sorted_evals,L,W):
     axs.format(
         xlim=(0, L), ylim=(0, -W), 
         )
+    
+def view_potential_shape(potential_type, n=2,Nc=1,a=6.75,Vo=0.306,closed_chain=False,L=100,W=100):
+    
+    lat = kwant.lattice.square(1)
+    syst = kwant.Builder()
+
+    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)
+                
+            elif potential_type == 'chain':
+                res = chain_n(pos, n=n, Nc=Nc, a=a, Vo=Vo, start_pos=[0,0], closed_chain=closed_chain)
+            if res == 0:
+                    return True
+            else:
+                return False
+
+    syst[lat.shape(shape_func, (0,0) )] = 1
+    syst[lat.neighbors(n=1)] = 1
+    kwant.plot(syst)
 
 

code/potentials_chain.py

@@ -5,6 +5,40 @@ import numpy as np
 import modules.potentials as potentials
 from matplotlib import pyplot as plt
 
+# # Patches
+
+# +
+hbar = 1
+m = 1
+a = 6.75 #0.36 nm (lattice constant) in a.u.
+
+
+lattice_spacing = 1
+L, W = 35, 35#sets the size of the discretized lattice (number points)
+well_pos = 2*a #sets the spacing between the wells
+
+
+chain = dict(
+    L=L,
+    W=W,
+    size=5,
+    start_pos = [0,0],
+    Nc = 0,
+    closed_chain = False,
+    a = a,
+    Vo = 0.306164,
+    potential_type = 'patch_n'
+)
+
+# +
+syst = potentials.potential_system(L,W,lattice_spacing,chain)
+sorted_evals, sorted_evecs = potentials.sorted_evs(syst, chain)
+
+x, y = np.linspace(0,L,L,endpoint=False), np.linspace(0, -W, W, endpoint=False)
+X, Y = np.meshgrid(x, y)
+potentials.plotter(X,Y,sorted_evecs,sorted_evals,L,W)
+# -
+
 # # Trimers
 
 # ## 2x2 trimer
@@ -28,6 +62,7 @@ chain = dict(
     Nc = 0,
     a = a,
     Vo = 0.306164,
+    closed_chain=False,
     potential_type = 'trimer_n'
 )
 
@@ -60,6 +95,7 @@ chain = dict(
     start_pos=[0,0],
     Nc = 0,
     Vo = 0.306164,
+    closed_chain=False,
     potential_type = 'trimer_n'
 )
 
@@ -94,6 +130,7 @@ chain = dict(
     start_pos=[0,0],
     Nc=0,
     Vo = 0.306164,
+    closed_chain=False,
     potential_type = 'trimer_n'
 )
 
@@ -130,6 +167,7 @@ chain = dict(
     start_pos=[0,0],
     Nc=2,
     Vo = 0.306164,
+    closed_chain=False,
     potential_type = 'chain_n'
 )
 
@@ -151,7 +189,7 @@ a = 6.75 #0.36 nm (lattice constant) in a.u.
 size = 3
 
 lattice_spacing = 1
-L, W = 150, 150 #sets the size of the discretized lattice (number points)
+L, W = 80, 80 #sets the size of the discretized lattice (number points)
 well_pos = 2*a #sets the spacing between the wells
 
 
@@ -161,8 +199,9 @@ chain = dict(
     a = a,
     size=size,
     start_pos=[0,0],
-    Nc=3,
+    Nc=2,
     Vo = 0.306164,
+    closed_chain=False,
     potential_type = 'chain_n'
 )
 
@@ -196,7 +235,8 @@ chain = dict(
     size=size,
     start_pos=[0,0],
     Nc=3,
-    Vo = 0.906164,
+    Vo = 0.306164,
+    closed_chain=True,
     potential_type = 'chain_n'
 )
 
@@ -208,5 +248,6 @@ x, y = np.linspace(0,L,L,endpoint=False), np.linspace(0, -W, W, endpoint=False)
 X, Y = np.meshgrid(x, y)
 potentials.plotter(X,Y,sorted_evecs,sorted_evals,L,W)
 # -
+potentials.view_potential_shape('chain',Nc=3,closed_chain=True,n=5,L=200,W=200)