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 Example usage:: 

35 

36 from tbkit.graphene import GrapheneLattice, GrapheneSystem 

37 

38 lat = GrapheneLattice() 

39 lat.hexagon_zigzag(n=4) 

40 sys = GrapheneSystem(lat) 

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

42 ''' 

43 

44 def __init__(self) -> None: 

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

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

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

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

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

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

51 

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

53 ''' 

54 Triangular flake with zigzag terminations. 

55 

56 :param: n. Int. Number of plaquettes along the edges.  

57 ''' 

58 error_handling.positive_int(n, 'n') 

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

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

61 

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

63 ''' 

64 Hexagonal flake with zigzag terminations. 

65 

66 :param: n. Int. Number of plaquettes along the edges.  

67 ''' 

68 error_handling.positive_int(n, 'n') 

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

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

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

72 

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

74 ''' 

75 Triangular flake with armchair terminations. 

76 

77 :param: n. Int. Number of plaquettes along the edges.  

78 ''' 

79 error_handling.positive_int(n, 'n') 

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

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

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

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

84 

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

86 ''' 

87 Hexagonal flake with armchair terminations. 

88 

89 :param: n. Int. Number of plaquettes along each edge.  

90 ''' 

91 error_handling.positive_int(n, 'n') 

92 nn = 3 * n - 2 

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

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

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

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

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

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

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

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

101 

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

103 ''' 

104 Squared flake. 

105 

106 :param: n. Int. Number of plaquettes along x.  

107 ''' 

108 error_handling.positive_int(n, 'n') 

109 n2 = int(1.5*DX*n) 

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

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

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

113 

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

115 ''' 

116 Circular flake. 

117 

118 :param: n. Int. Number of plaquettes along the diameter.  

119 ''' 

120 error_handling.positive_int(n, 'n') 

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

122 self.sites = len(self.coor) 

123 self.center() 

124 if n % 2 == 0: 

125 self.shift_x(shift=-DX) 

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

127 self.remove_dangling() 

128 

129class GrapheneSystem(System): 

130 r''' 

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

132 

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

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

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

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

137 (*get_butterfly*). 

138 

139 :param lat: Lattice, usually a :class:`GrapheneLattice` flake. 

140 

141 Example usage:: 

142 

143 from tbkit.graphene import GrapheneLattice, GrapheneSystem 

144 

145 lat = GrapheneLattice() 

146 lat.hexagon_zigzag(n=10) 

147 lat.center() 

148 sys = GrapheneSystem(lat) 

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

150 sys.get_ham() 

151 ''' 

152 

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

154 System.__init__(self, lat) 

155 

156 def _strain_projection( 

157 self, 

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

159 r''' 

160 Private method. 

161 

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

163 

164 .. math:: 

165 

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

167 

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

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

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

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

172 

173 :returns: 

174 * **i**, **j** -- Integer ndarrays. Bond end points, oriented as in 

175 *System.get_bonds*. 

176 * **ang** -- Real ndarray. Bond angles (from *i* to *j*), in degrees, 

177 in [0, 180). 

178 * **s** -- Real ndarray. The projections above. 

179 ''' 

180 self.get_distances() 

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

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

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

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

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

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

187 return i, j, ang, s 

188 

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

190 r''' 

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

192 

193 .. math:: 

194 

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

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

197 

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

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

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

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

202 before calling this. 

203 

204 :param t: Hopping value without strain. 

205 :param beta: Strength of the strain. See *get_beta_lims* for the 

206 range that keeps every hopping positive. 

207 ''' 

208 error_handling.number(t, 't') 

209 error_handling.real_number(beta, 'beta') 

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

211 self.clear_hopping() 

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

213 self.hop['n'] = 1 

214 self.hop['i'] = i 

215 self.hop['j'] = j 

216 self.hop['ang'] = ang 

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

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

219 

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

221 ''' 

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

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

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

225 

226 :param t: Unstrained hopping value. 

227 :param N: Positive integer. Number of strain values between the 

228 minimal and maximal strains given by *get_beta_lims*. 

229 ''' 

230 error_handling.number(t, 't') 

231 error_handling.positive_int(N, 'N') 

232 beta_lims = self.get_beta_lims() 

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

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

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

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

237 self.get_ham() 

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

239 

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

241 r''' 

242 Get the extremal strain values keeping every hopping positive. 

243 

244 Under *set_hop_linear_strain* a bond's amplitude is 

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

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

247 

248 :returns: 

249 * **beta_lims** -- Real ndarray of length 2, ``[beta_min, beta_max]`` 

250 (ascending). A bound is infinite if the corresponding 

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

252 ''' 

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

254 eps = 1e-6 

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

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

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

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

259 

260 

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

262grapheneLat = GrapheneLattice 

263grapheneSys = GrapheneSystem