Coverage for tbkit/graphene.py: 100%

1from __future__ import annotations 

2 

3from tbkit.lattice import * 

4from tbkit.plot import * 

5from tbkit.system import * 

6from math import sqrt 

7 

8PI = np.pi 

9ATOL = 1e-3 

10DX = 0.5 * sqrt(3) 

11DY = 0.5 

12 

13 

14################################# 

15# CLASS GRAPHENE 

16################################# 

17 

18 

19class GrapheneLattice(Lattice): 

20 r''' 

21 Honeycomb lattice of graphene, with flake shapes. 

22 

23 Two sublattices, 'a' at :math:`(0, 0)` and 'b' at 

24 :math:`(\sqrt3/2, 1/2)`, and primitive vectors 

25 :math:`\mathbf{a}_1 = (\sqrt3, 0)` and 

26 :math:`\mathbf{a}_2 = (\sqrt3/2, 3/2)`, so the nearest-neighbour 

27 distance is 1. It is a :class:`Lattice` with methods that cut flakes 

28 of the usual shapes: *triangle_zigzag*, *hexagon_zigzag*, 

29 *triangle_armchair*, *hexagon_armchair*, *square* and *circle*. 

30 

31 *butterfly* and *betas* are empty here; *GrapheneSystem.get_butterfly* 

32 fills the ones of the system. 

33 

34 Examples 

35 -------- 

36 Example usage:: 

37 

38 from tbkit.graphene import GrapheneLattice, GrapheneSystem 

39 

40 lat = GrapheneLattice() 

41 lat.hexagon_zigzag(n=4) 

42 sys = GrapheneSystem(lat) 

43 sys.set_hopping([{'n': 1, 't': 1.}]) 

44 ''' 

45 

46 def __init__(self) -> None: 

47 unit_cell = [{'tag': 'a', 'r0': (0, 0)}, 

48 {'tag': 'b', 'r0': (DX, DY)}] 

49 prim_vec = [(2*DX, 0.), (DX, 1.5)] 

50 Lattice.__init__(self, unit_cell=unit_cell, prim_vec=prim_vec) 

51 self.butterfly = np.array([]) 

52 self.betas = np.array([]) 

53 

54 def triangle_zigzag(self, n: int) -> None: 

55 ''' 

56 Triangular flake with zigzag terminations. 

57 

58 Parameters 

59 ---------- 

60 n : Int 

61 Number of plaquettes along the edges. 

62 ''' 

63 error_handling.positive_int(n, 'n') 

64 self.get_lattice(n1=n+2, n2=n+2) 

65 self.boundary_line(cx=-sqrt(3), cy=-1, co=-3*n+2) 

66 

67 def hexagon_zigzag(self, n: int) -> None: 

68 ''' 

69 Hexagonal flake with zigzag terminations. 

70 

71 Parameters 

72 ---------- 

73 n : Int 

74 Number of plaquettes along the edges. 

75 ''' 

76 error_handling.positive_int(n, 'n') 

77 self.get_lattice(n1=2*n, n2=2*n) 

78 self.boundary_line(cx=sqrt(3), cy=1, co=3*(n-1)) 

79 self.boundary_line(cx=-sqrt(3), cy=-1, co=-9*n+2.5) 

80 

81 def triangle_armchair(self, n: int) -> None: 

82 ''' 

83 Triangular flake with armchair terminations. 

84 

85 Parameters 

86 ---------- 

87 n : Int 

88 Number of plaquettes along the edges. 

89 ''' 

90 error_handling.positive_int(n, 'n') 

91 self.get_lattice(n1=2*n, n2=2*n) 

92 self.boundary_line(cx=1, cy=0, co=sqrt(3)/2*(2*n-1)-0.1) 

93 self.boundary_line(cx=-1/sqrt(3), cy=1, co=-n-0.1) 

94 self.boundary_line(cx=-1/sqrt(3), cy=-1, co=-4*n+0.1) 

95 

96 def hexagon_armchair(self, n: int) -> None: 

97 ''' 

98 Hexagonal flake with armchair terminations. 

99 

100 Parameters 

101 ---------- 

102 n : Int 

103 Number of plaquettes along each edge. 

104 ''' 

105 error_handling.positive_int(n, 'n') 

106 nn = 3 * n - 2 

107 self.get_lattice(n1=2*nn, n2=2*nn) 

108 self.boundary_line(cx=1, cy=0, co=sqrt(3)/2* (2*nn-1)-.1) 

109 self.boundary_line(cx=-1/sqrt(3), cy=1, co=-nn -.1) 

110 self.boundary_line(cx=-1/sqrt(3), cy=-1, co=-4*nn +1-.1) 

111 self.coor['x'] -= self.coor['x'].min() 

112 self.boundary_line(cx=-1, cy=0, co=-DX * 2*nn-0.1) 

113 self.boundary_line(cx=1/sqrt(3), cy=1, co=3*n -3.) 

114 self.boundary_line(cx=1/sqrt(3), cy=-1, co=-6*n+4) 

115 

116 def square(self, n: int) -> None: 

117 ''' 

118 Squared flake. 

119 

120 Parameters 

121 ---------- 

122 n : Int 

123 Number of plaquettes along x. 

124 ''' 

125 error_handling.positive_int(n, 'n') 

126 n2 = int(1.5*DX*n) 

127 self.get_lattice(n1=2*n, n2=n2) 

128 self.boundary_line(cx=1, cy=0, co=DX*(2*n-2)+0.1) 

129 self.boundary_line(cx=-1, cy=0, co=-DX*(4*n)+0.5) 

130 

131 def circle(self, n: int) -> None: 

132 ''' 

133 Circular flake. 

134 

135 Parameters 

136 ---------- 

137 n : Int 

138 Number of plaquettes along the diameter. 

139 ''' 

140 error_handling.positive_int(n, 'n') 

141 self.get_lattice(n1=2*n, n2=2*n) 

142 self.sites = len(self.coor) 

143 self.center() 

144 if n % 2 == 0: 

145 self.shift_x(shift=-DX) 

146 self.ellipse_in(rx=DX*(n+1), ry=DX*(n+1), x0=0., y0=0.) 

147 self.remove_dangling() 

148 

149class GrapheneSystem(System): 

150 r''' 

151 Graphene :class:`System`, with triaxial strain. 

152 

153 Adds to *System* the nearest-neighbour hoppings of a linear triaxial 

154 strain (*set_hop_linear_strain*), which acts as a uniform 

155 pseudo-magnetic field, the range of strains that keeps every hopping 

156 positive (*get_beta_lims*), and the spectrum as a function of the strain 

157 (*get_butterfly*). 

158 

159 Parameters 

160 ---------- 

161 lat : Lattice, usually a :class:`GrapheneLattice` flake 

162 

163 Examples 

164 -------- 

165 Example usage:: 

166 

167 from tbkit.graphene import GrapheneLattice, GrapheneSystem 

168 

169 lat = GrapheneLattice() 

170 lat.hexagon_zigzag(n=10) 

171 lat.center() 

172 sys = GrapheneSystem(lat) 

173 sys.set_hop_linear_strain(t=1., beta=0.1) 

174 sys.get_ham() 

175 ''' 

176 

177 def __init__(self, lat: Lattice) -> None: 

178 System.__init__(self, lat) 

179 

180 def _strain_projection( 

181 self, 

182 ) -> tuple[NDArray[np.int64], NDArray[np.int64], NDArray[np.float64], NDArray[np.float64]]: 

183 r''' 

184 Private method. 

185 

186 Get, for every nearest-neighbor bond, the projection 

187 

188 .. math:: 

189 

190 s_{ij} = \hat{\boldsymbol\delta}_{ij}\cdot\mathbf{r}_{ij} 

191 

192 of its midpoint :math:`\mathbf{r}_{ij}` on its direction 

193 :math:`\hat{\boldsymbol\delta}_{ij}`, oriented from its 'b' site to 

194 its 'a' site (the same orientation for all three bond families), so 

195 that the strained hopping is :math:`t_{ij} = t(1 + \tfrac14\beta s_{ij})`. 

196 

197 Returns 

198 ------- 

199 i, j : Integer ndarrays 

200 Bond end points, oriented as in 

201 *System.get_bonds*. 

202 ang : Real ndarray 

203 Bond angles (from *i* to *j*), in degrees, 

204 in [0, 180). 

205 s : Real ndarray 

206 The projections above. 

207 ''' 

208 self.get_distances() 

209 i, j, ang = self.get_bonds(1) 

210 # direction from the 'b' end to the 'a' end of each bond 

211 ang_ba = np.where(self.lat.coor['tag'][i] == 'a', ang - 180., ang) 

212 x_center = .5 * (self.lat.coor['x'][i] + self.lat.coor['x'][j]) 

213 y_center = .5 * (self.lat.coor['y'][i] + self.lat.coor['y'][j]) 

214 s = (np.cos(PI / 180 * ang_ba) * x_center + np.sin(PI / 180 * ang_ba) * y_center) 

215 return i, j, ang, s 

216 

217 def set_hop_linear_strain(self, t: complex, beta: float) -> None: 

218 r''' 

219 Set nearest-neighbor hoppings according to a linear triaxial strain: 

220 

221 .. math:: 

222 

223 t_{ij} = t\left(1 + \tfrac14\beta\, 

224 \hat{\boldsymbol\delta}_{ij}\cdot\mathbf{r}_{ij}\right) 

225 

226 with :math:`\hat{\boldsymbol\delta}_{ij}` the bond direction (from its 

227 'b' site to its 'a' site) and 

228 :math:`\mathbf{r}_{ij}` its midpoint. The strain is measured from the 

229 coordinate origin, so centre the flake on it (see *lattice.center*) 

230 before calling this. 

231 

232 Parameters 

233 ---------- 

234 t 

235 Hopping value without strain. 

236 beta 

237 Strength of the strain. See *get_beta_lims* for the 

238 range that keeps every hopping positive. 

239 ''' 

240 error_handling.number(t, 't') 

241 error_handling.real_number(beta, 'beta') 

242 i, j, ang, s = self._strain_projection() 

243 self.clear_hopping() 

244 self.hop = np.zeros(len(i), dtype=HOP_DTYPE) 

245 self.hop['n'] = 1 

246 self.hop['i'] = i 

247 self.hop['j'] = j 

248 self.hop['ang'] = ang 

249 self.hop['tag'] = npc.add(self.lat.coor['tag'][i], self.lat.coor['tag'][j]) 

250 self.hop['t'] = t * (1. + 0.25 * beta * s) 

251 

252 def get_butterfly(self, t: complex, N: int) -> None: 

253 ''' 

254 Get the spectrum as a function of the strain *beta* (see 

255 *set_hop_linear_strain*), stored in *sys.butterfly*, shape 

256 (N, sites), with the strain values in *sys.betas*. 

257 

258 Parameters 

259 ---------- 

260 t 

261 Unstrained hopping value. 

262 N : Positive integer 

263 Number of strain values between the 

264 minimal and maximal strains given by *get_beta_lims*. 

265 ''' 

266 error_handling.number(t, 't') 

267 error_handling.positive_int(N, 'N') 

268 beta_lims = self.get_beta_lims() 

269 self.betas = np.linspace(beta_lims[0], beta_lims[1], N) 

270 self.butterfly = np.zeros((N, self.lat.sites)) 

271 for i, beta in enumerate(self.betas): 

272 self.set_hop_linear_strain(t=t, beta=beta) 

273 self.get_ham() 

274 self.butterfly[i] = LA.eigvalsh(self.ham.toarray()) 

275 

276 def get_beta_lims(self) -> NDArray[np.float64]: 

277 r''' 

278 Get the extremal strain values keeping every hopping positive. 

279 

280 Under *set_hop_linear_strain* a bond's amplitude is 

281 :math:`t(1+\tfrac14\beta s_{ij})`, so it stays positive for every bond 

282 iff :math:`-4/\max_{ij} s_{ij} < \beta < -4/\min_{ij} s_{ij}`. 

283 

284 Returns 

285 ------- 

286 beta_lims : Real ndarray of length 2, ``[beta_min, beta_max]`` (ascending) 

287 A bound is infinite if the corresponding 

288 :math:`s_{ij}` never takes that sign. 

289 ''' 

290 _, _, _, s = self._strain_projection() 

291 eps = 1e-6 

292 s_min, s_max = s.min(), s.max() 

293 beta_min = -4. / s_max + eps if s_max > 0 else -np.inf 

294 beta_max = -4. / s_min - eps if s_min < 0 else np.inf 

295 return np.array([beta_min, beta_max]) 

296 

297 

298# Backward-compatible camelCase aliases (pre-0.2 API). 

299grapheneLat = GrapheneLattice 

300grapheneSys = GrapheneSystem