-
Notifications
You must be signed in to change notification settings - Fork 1
/
HyperBdG
201 lines (158 loc) · 6.53 KB
/
HyperBdG
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
import numpy as np
from scipy.linalg import eigh
import matplotlib.pyplot as plt
from matplotlib import cm, rc, ticker
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
from scipy.sparse import diags
from numpy import random
import pickle
import time
import matplotlib.colors as mcolors
import networkx as nx
import math
"Class for initial one-particle lattice and corresponding hamiltonian"
class HyperLattice:
def __init__(self,p,q,l,hopping):
self.p=p
self.q=q
self.l=l
self.hopping=hopping
self.sites=np.array([])
self.hamiltonian=[]
self.create_hyperbolic_lattice()
def gamma(self, z: np.complex128) -> np.complex128:
sigma = math.sqrt( (math.cos(2*math.pi / p) + math.cos(2*math.pi / q)) / (1 + math.cos(2*math.pi / q)) )
z = np.complex128(z)
result = (z + sigma) / (sigma * z + 1)
return result
def rot(self, z, n):
z = np.complex128(z)
result = (math.cos(2 * math.pi * n / p) + math.sin(2 * math.pi * n / p) * 1j) * z
return result
def trans(self, z, n):
result = self.rot(z,-n)
result = self.gamma(result)
result = self.rot(result, n)
return result
# It is useful to define the distance function
def dist(self, x, y):
return math.acosh(1 + 2*abs(x-y)**2 / (1 - abs(x)**2) / (1-abs(y)**2))
def create_hyperbolic_lattice(self):
r0 = math.sqrt( math.cos(math.pi*(1/self.p +1/self.q) ) / math.cos(math.pi*(1/self.p - 1/self.q) ) )
#Unit cell points
for n in range(p):
self.sites = np.append(self.sites, r0*math.cos(math.pi*(1+2*n)/self.p) + r0*math.sin(math.pi*(1+2*n)/self.p)*1j )
# Next, we generate new cells by applying translation generators (and their inverses).
i=1
#print(sites.size)
while i < self.l:
i=i+1
N=self.sites.size
for k in range(N):
for n in range(N):
self.sites = np.append(self.sites, self.trans(self.sites[k], n)) # we apply generators to each side...
self.sites = np.unique(self.sites) # and through away the repeated ones.
#self.sites = np.unique(self.sites)
#centers = np.array([])
#i=1
#while i < l:
# i=i+1
# for n in range(p):
# centers = np.append(centers, trans(0, n)) # we apply generators to each side...
# Let us check again that no repeated sites are generated, and if they are, we through them away.
i=0
while i < self.sites.size:
ind_to_del = np.array([], dtype=int)
for k in range(i+1, self.sites.size):
if self.dist(self.sites[i], self.sites[k]) < 0.01:
ind_to_del=np.append(ind_to_del, k)
self.sites = np.delete(self.sites, ind_to_del)
i=i+1
# Having generated the lattice, we can now build the adjacncy matrix.
print('We have ', self.sites.size,' sites in total.')
N=self.sites.size
C = self.dist(r0, 0)
B = math.asinh( math.sin( math.pi / p )*math.sinh(C))
self.hamiltonian=np.zeros((N,N))
for i in range(N):
for k in range(N):
if self.dist(self.sites[i], self.sites[k]) < 2*B+0.001:
if self.dist(self.sites[i], self.sites[k]) > 2*B-0.001:
self.hamiltonian[i, k] = -self.hopping
print(self.hamiltonian)
"Class for BdG hamiltonians and all corresponding functions"
class HyperBdG():
def __init__(self, hyperlattice, V,T,mu,Delta=[]):
self.lattice_sample=hyperlattice
self.lattice_H=hyperlattice.hamiltonian
self.N=len(hyperlattice.sites)
self.hopping=hyperlattice.hopping
self.V=V
self.T=T
self.mu=mu
self.BdG_H=[]
if len(Delta)==0:
self.Delta=np.zeros(self.N)
self.initial_Delta=False #it is necessary to obtain non-trivial solution of BdG equation
print("no initial Delta for the BdG Hamiltonian")
else:
self.Delta=Delta
self.initial_Delta=True
#construct trivial hamiltonian
self.construct_hamiltonian()
spectra, vectors = eigh(self.BdG_H)
self.spectra=spectra
self.vectors=vectors
#Fermi function
def F(self, E):
return 1/(np.exp((E)/self.T)+1)
def construct_hamiltonian(self):
H_Delta = diags([self.Delta], [0], shape=(self.N, self.N)).toarray()
self.BdG_H = np.block([[self.lattice_H - self.mu*np.eye(self.N), H_Delta], [H_Delta, -self.lattice_H + self.mu*np.eye(self.N)]])
def charge_density(self):
fermi_dist=self.F(self.spectra[self.N:])
v=self.vectors[self.N:,self.N:]
u=self.vectors[:self.N,self.N:]
n=2*np.einsum(u,[0,1],np.conj(u), [0,1],fermi_dist,[1],[0])+2*np.einsum(v,[0,1],np.conj(v), [0,1],np.ones(self.N)-fermi_dist,[1],[0])
return n
def BdG_cycle(self):
print("charge density, n", np.mean(self.charge_density()))
print("BdG cycle T=", self.T)
step=0
if self.initial_Delta==False:
self.Delta=0.5*np.ones(self.N)+0.1*np.random.rand(self.N)
self.construct_hamiltonian()
spectra, vectors = eigh(self.BdG_H)
self.spectra=spectra
self.vectors=vectors
while True:
F_weight=np.ones(self.N)-2*self.F(self.spectra[self.N:])
vectors_up=self.V * np.conj(self.vectors[self.N:,self.N:])
Delta_next= np.einsum(vectors_up, [0,1], self.vectors[:self.N,self.N:], [0,1], F_weight,[1],[0])
error=np.max(np.abs((self.Delta-Delta_next)))
self.Delta=Delta_next
print("step", step, "error", error, "Delta_max", np.max(np.abs(self.Delta)))
step += 1
if error<10**(-6):
break
self.construct_hamiltonian()
spectra, vectors = eigh(self.BdG_H)
self.spectra=spectra
self.vectors=vectors
"Sample run input"
p=7
q=3
l=3
t=1
V=5
T=0.02
mu=0
hypersample=HyperLattice(p,q,l,t)
#print(hypersample.hamiltonian)
#print(hypersample.sites)
BdGhypersample=HyperBdG(hypersample,V,T,mu)
BdGhypersample.BdG_cycle()
print(BdGhypersample.Delta)