Implement cross-coupling between lattices of different helicities via the coupling patches

code/modules/potentials.py

@@ -275,6 +275,89 @@ def onsite(site, Vo, potential_type, a, size, start_pos, Nc, closed_chain, inver
     elif potential_type == 'ssh_n':
         return 4*t + ssh_n(site, size, a, Vo, start_pos,trivial,dim=dim)
     
+    
+def hopping(site1,site2, potential_type, size, Nc, a, Vo, dim, closed_chain, inversion_symmetric, start_pos, cross_hop_x,cross_hop_y):
+    
+    x1,y1 = site1.pos
+    x2,y2 = site2.pos
+    
+    if np.abs(x1-x2) > 1 or np.abs(y1-y2):
+        return 0
+    
+    if potential_type == 'trimer_n':
+        Nc_ = 1
+    elif potential_type == 'chain_n':
+        Nc_ = Nc
+    
+    #connecting patch size
+    n = size
+    n_ = n-1-n%2
+    
+    if inversion_symmetric:
+        n__ = False
+    elif not inversion_symmetric and a > 1:
+        n__ = n_*a//2
+        
+    def in_patch(site):
+    
+        for Uc in range(Nc_):
+            if Uc != 0:
+                trimer_start_pos = [start_pos[0]+Uc*a*(n+n_), start_pos[1]+Uc*a*(n+n_)]
+            else:
+                trimer_start_pos = start_pos
+
+
+            if n%2 == 0:
+                patch_1_2 = patch_n(site,a, Vo, n_, [trimer_start_pos[0]+n*a, trimer_start_pos[1]],rectangle=n__,dim=dim)
+                if patch_1_2 == 0:
+                    return 'x'
+                patch_1_3 = patch_n(site,a, Vo, n_, [trimer_start_pos[0], trimer_start_pos[1]+n*a],dim=dim)
+                if patch_1_3 == 0:
+                    return 'y'
+                if Uc!=0:
+                    extra_connector_1_2 = patch_n(site,a, Vo, n_, [trimer_start_pos[0]+1*a, trimer_start_pos[1]-n_*a],dim=dim)
+                    if extra_connector_1_2 == 0:
+                        return 'y'
+                    extra_connector_1_3 = patch_n(site,a, Vo, n_, [trimer_start_pos[0]-n_*a, trimer_start_pos[1]+1*a],dim=dim)
+                    if extra_connector_1_3 == 0:
+                        return 'x'
+
+            else:
+                patch_1_2 = patch_n(site,a, Vo, n_, [trimer_start_pos[0]+n*a, trimer_start_pos[1]+1*a],rectangle=n__,dim=dim)
+                if patch_1_2 == 0:
+                    return 'x'
+                patch_1_3 = patch_n(site,a, Vo, n_, [trimer_start_pos[0]+1*a, trimer_start_pos[1]+n*a],dim=dim)
+                if patch_1_3 == 0:
+                    return 'y'
+                if Uc != 0:
+                    extra_connector_1_2 = patch_n(site,a, Vo, n_, [trimer_start_pos[0]+1*a, trimer_start_pos[1]-n_*a],dim=dim)
+                    if extra_connector_1_2 == 0:
+                        return 'y'
+                    extra_connector_1_3 = patch_n(site,a, Vo, n_, [trimer_start_pos[0]-n_*a, trimer_start_pos[1]+1*a],dim=dim)
+                    if extra_connector_1_3 == 0:
+                        return 'x'
+                    
+        return 
+    
+    def hopping_per_direction(direction):
+        if direction == 'x':
+            return cross_hop_x
+        elif direction == 'y':
+            return cross_hop_y
+        else:
+            return 0
+                    
+    in_patch_1 = in_patch(site1)
+    if in_patch_1:
+        in_patch_2 = in_patch(site2)
+        if in_patch_2:
+            if in_patch_1 != in_patch_2:
+                raise runtimeError('different connecting patches!')
+            return hopping_per_direction(in_patch_1)
+        else:
+            return 0
+    else:
+        return 0
 
 
 # -
@@ -291,28 +374,45 @@ def potential_system(L,W,params,dim,cross_hop_x, cross_hop_y):
     if dim == 1:
         lat = kwant.lattice.square(1,norbs=1) #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(n=1)] = np.eye(dim)*t
+        syst[lat.neighbors(n=1)] = t
         
     elif dim == 2:
         lat_1 = kwant.lattice.square(1, name='1', norbs=1)
         lat_2 = kwant.lattice.square(1, name='2',norbs=1)
-        syst[(lat_1(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_2(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_1.neighbors(n=1)] = t
-#         syst[lat_2.neighbors(n=1)] = t
         
-        syst[kwant.builder.HoppingKind((1, 0), lat_1, lat_1)] = t
-        syst[kwant.builder.HoppingKind((0, 1), lat_1, lat_1)] = t
         
+        for i in range(L):
+            for j in np.linspace(0,-W,W,endpoint=False):
+                syst[lat_1(i, j)] = onsite
+                syst[lat_2(i, j)] = onsite
+                
+        hoppings = (((0, 1), lat_1, lat_2), ((1, 0), lat_1, lat_2))
+        syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = hopping
         
-        syst[kwant.builder.HoppingKind((1, 0), lat_2, lat_2)] = t
-        syst[kwant.builder.HoppingKind((0, 1), lat_2, lat_2)] = t
         
-
-        syst[kwant.builder.HoppingKind((1, 0), lat_1, lat_2)] = cross_hop_x
-        syst[kwant.builder.HoppingKind((1, 0), lat_2, lat_1)] = cross_hop_x
-        syst[kwant.builder.HoppingKind((0, 1), lat_1, lat_2)] = cross_hop_y
-        syst[kwant.builder.HoppingKind((0, 1), lat_2, lat_1)] = cross_hop_y
+#         for i in range(L):
+#             for j in np.linspace(0,-W,W,endpoint=False):
+                
+#                 i_lt_0 = i>0
+#                 i_st_L = i < L-1
+#                 j_st_0 = j < 0
+#                 j_lt_W = j > -W+1
+               
+                
+#                 if i_lt_0:
+#                     syst[lat_1(i,j),lat_2(i-1,j)] = hopping
+#                 if i_st_L:
+#                     syst[lat_1(i,j),lat_2(i+1,j)] = hopping
+#                 if j_st_0:
+#                     syst[lat_1(i,j),lat_2(i,j+1)] = hopping
+#                 if j_lt_W:
+#                     syst[lat_1(i,j),lat_2(i,j-1)] = hopping
+                    
+                    
+                    
+        syst[lat_1.neighbors(n=1)] = t
+        syst[lat_2.neighbors(n=1)] = t
+                    
 
     syst = syst.finalized() #hamiltonian is built
 

code/potentials_chain.py

@@ -184,7 +184,7 @@ a = 6.75 #0.36 nm (lattice constant) in a.u.
 lattice_spacing = 1
 
 size = 3
-Nc = 1
+Nc = 3
 
 
 size_ = size-1-size%2
@@ -200,16 +200,18 @@ chain = dict(
     size=size,
     start_pos=[0,0],
     Nc=Nc,
-    Vo = 0.0506164,
+    Vo = 0.506164,
     dim=2,
     closed_chain=False,
     inversion_symmetric=True,
     potential_type = 'chain_n',
+    cross_hop_x = 0.,
+    cross_hop_y = -0.,
     trivial=False
 )
 
 # +
-syst = potentials.potential_system(L,W,chain,dim=chain['dim'],cross_hop_x=0.5,cross_hop_y=-0.5)
+syst = potentials.potential_system(L,W,chain, dim=chain['dim'], cross_hop_x=chain['cross_hop_x'],cross_hop_y=chain['cross_hop_y'])
 # syst = potentials.view_potential_shape(L=chain['L']
 #                             ,W=chain['W'],
 #                             n=chain['size'],
@@ -223,15 +225,15 @@ syst = potentials.potential_system(L,W,chain,dim=chain['dim'],cross_hop_x=0.5,cr
 #                             cross_coupling_y=0,
 #                             semi_infinite=False)
 
-sorted_evals, sorted_evecs = potentials.sorted_evs(syst, chain,k=100)
+sorted_evals, sorted_evecs = potentials.sorted_evs(syst, chain,k=30)
 
 # -
-sorted_evals
-
 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,rounding=2)
-potentials.plotter(X,Y,sorted_evecs[len(sorted_evecs)//2:,:],sorted_evals,L,W,rounding=3)
+if chain['dim'] == 1:
+    potentials.plotter(X,Y,sorted_evecs,sorted_evals,L,W,rounding=2)
+elif chain['dim'] == 2:
+    potentials.plotter(X,Y,sorted_evecs[len(sorted_evecs)//2:,:],sorted_evals,L,W,rounding=2)
 
 # # SSH chain
 
@@ -277,124 +279,3 @@ 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, aspect = 2)
-# -
-
-# # Complex system
-
-# +
-hbar = 1
-m = 1
-a = 6.75 #0.36 nm (lattice constant) in a.u.
-lattice_spacing = 1
-
-size = 3
-Nc = 2
-
-
-size_ = size-1-size%2
-L = int( (size+size_)*a*Nc + size*a )
-W = L#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,
-    a = a,
-    size=size,
-    start_pos=[0,0],
-    Nc=Nc,
-    Vo = 0.506164,
-    phi=np.pi/2,
-    closed_chain=False,
-    inversion_symmetric=False,
-    potential_type = 'chain_n',
-    trivial=False,
-)
-
-# +
-syst = potentials.complex_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)
-# -
-
-def sorted_evs(syst, potential):
-    '''
-    :param kwant.Builder syst: kwant system of substrate
-    '''
-    
-    ham_mat = syst.hamiltonian_submatrix(params=potential, sparse=True)
-    evals, evecs = sla.eigsh(ham_mat.tocsc(), k=2, sigma=0, return_eigenvectors=True)
-    
-    evals = evals*27.211324570273+4.7
-    inds = evals.argsort()
-    sorted_evals = evals[inds]
-    sorted_evecs = evecs[:,inds]
-    return evals, evecs, sorted_evals, sorted_evecs
-
-
-ham_mat = syst.hamiltonian_submatrix(params=chain, sparse=True)
-
-sla.eigsh(ham_mat.toarray(), k=2, sigma=0, return_eigenvectors=True)
-
-ham_mat.shape
-
-plt.spy(ham_mat)
-
-test = np.nonzero(np.round(sevecs,2))[0]
-
-test
-
-L*6
-
-np.abs(sevecs[test].reshape((L*6, W)).T)**2
-
-sevecs[L*W:L*W+5,0]
-
-
-def plotter_6_sublattices(X,Y,sorted_evecs,sorted_evals,L,W, aspect = None):
-    
-    column_no = 8
-    row_no = len(sorted_evals)//8+1
-    
-    fig, axs = plot.subplots(ncols=column_no, nrows=row_no)
-
-    props = dict(N=100, rasterized=True, symmetric=True, cmap='dusk_r')
-
-    j=0
-    for row in range(row_no):
-        for i in range(column_no):
-            try:
-                b = sorted_evals[j+1]
-            except:
-                continue
-                
-            for divide_in_6 in range(6):
-                evals
-                
-            k = 0
-            while round(sorted_evals[j+k],2) == round(sorted_evals[j+1+k],2):
-                k += 1
-                
-            colormesh_ = np.zeros((L*6,W)).T
-            for k_i in range(k):
-                colormesh_+= np.abs(sorted_evecs[:,j+k_i].reshape((L*6, W)).T)**2
-
-                
-            axs[row,i].pcolormesh(X, Y, colormesh_, **props)
-            j=j+k
-            
-            axs.format(
-                xlim=(0, L), ylim=(0, -W), 
-                )
-            if aspect:
-                axs.set_aspect(aspect)
-
-plotter_6_sublattices(None,None,sorted_evecs,sorted_evals,L,W)
-
-74*74*6
-
-