1"""Bravais.py - class for generating Bravais lattices etc. 

2 

3This is a base class for numerous classes setting up pieces of crystal. 

4""" 

5 

6import math 

7from typing import Optional, Sequence 

8 

9import numpy as np 

10 

11from ase.atoms import Atoms 

12import ase.data 

13 

14 

15class Bravais: 

16 """Bravais lattice factory. 

17 

18 This is a base class for the objects producing various lattices 

19 (SC, FCC, ...). 

20 """ 

21 

22 # The following methods are NOT defined here, but must be defined 

23 # in classes inhering from Bravais: 

24 # get_lattice_constant 

25 # make_crystal_basis 

26 # The following class attributes are NOT defined here, but must be defined 

27 # in classes inhering from Bravais: 

28 # int_basis 

29 # inverse_basis 

30 

31 other = {0: (1, 2), 1: (2, 0), 2: (0, 1)} 

32 

33 # For Bravais lattices with a basis, set the basis here. Leave as 

34 # None if no basis is present. 

35 bravais_basis: Optional[Sequence[Sequence[float]]] = None 

36 

37 # If more than one type of element appear in the crystal, give the 

38 # order here. For example, if two elements appear in a 3:1 ratio, 

39 # bravais_basis could contain four vectors, and element_basis 

40 # could be (0,0,1,0) - the third atom in the basis is different 

41 # from the other three. Leave as None if all atoms are of the 

42 # same type. 

43 element_basis: Optional[Sequence[int]] = None 

44 

45 # How small numbers should be considered zero in the unit cell? 

46 chop_tolerance = 1e-10 

47 

48 def __call__(self, symbol, 

49 directions=(None, None, None), miller=(None, None, None), 

50 size=(1, 1, 1), latticeconstant=None, 

51 pbc=True, align=True, debug=0): 

52 "Create a lattice." 

53 self.size = size 

54 self.pbc = pbc 

55 self.debug = debug 

56 self.process_element(symbol) 

57 self.find_directions(directions, miller) 

58 if self.debug: 

59 self.print_directions_and_miller() 

60 self.convert_to_natural_basis() 

61 if self.debug >= 2: 

62 self.print_directions_and_miller(" (natural basis)") 

63 if latticeconstant is None: 

64 if self.element_basis is None: 

65 self.latticeconstant = self.get_lattice_constant() 

66 else: 

67 raise ValueError( 

68 "A lattice constant must be specified for a compound") 

69 else: 

70 self.latticeconstant = latticeconstant 

71 if self.debug: 

72 print( 

73 "Expected number of atoms in unit cell:", 

74 self.calc_num_atoms()) 

75 if self.debug >= 2: 

76 print("Bravais lattice basis:", self.bravais_basis) 

77 if self.bravais_basis is not None: 

78 print(" ... in natural basis:", self.natural_bravais_basis) 

79 self.make_crystal_basis() 

80 self.make_unit_cell() 

81 if align: 

82 self.align() 

83 return self.make_list_of_atoms() 

84 

85 def align(self): 

86 "Align the first axis along x-axis and the second in the x-y plane." 

87 degree = 180 / np.pi 

88 if self.debug >= 2: 

89 print("Basis before alignment:") 

90 print(self.basis) 

91 if self.basis[0][0]**2 + \ 

92 self.basis[0][2]**2 < 0.01 * self.basis[0][1]**2: 

93 # First basis vector along y axis - rotate 90 deg along z 

94 t = np.array([[0, -1, 0], 

95 [1, 0, 0], 

96 [0, 0, 1]], float) 

97 self.basis = np.dot(self.basis, t) 

98 transf = t 

99 if self.debug >= 2: 

100 print("Rotating -90 degrees around z axis for numerical " 

101 "stability.") 

102 print(self.basis) 

103 else: 

104 transf = np.identity(3, float) 

105 assert abs(np.linalg.det(transf) - 1) < 1e-6 

106 # Rotate first basis vector into xy plane 

107 theta = math.atan2(self.basis[0, 2], self.basis[0, 0]) 

108 t = np.array([[np.cos(theta), 0, -np.sin(theta)], 

109 [0, 1, 0], 

110 [np.sin(theta), 0, np.cos(theta)]]) 

111 self.basis = np.dot(self.basis, t) 

112 transf = np.dot(transf, t) 

113 if self.debug >= 2: 

114 print("Rotating %f degrees around y axis." % (-theta * degree,)) 

115 print(self.basis) 

116 assert abs(np.linalg.det(transf) - 1) < 1e-6 

117 # Rotate first basis vector to point along x axis 

118 theta = math.atan2(self.basis[0, 1], self.basis[0, 0]) 

119 t = np.array([[np.cos(theta), -np.sin(theta), 0], 

120 [np.sin(theta), np.cos(theta), 0], 

121 [0, 0, 1]]) 

122 self.basis = np.dot(self.basis, t) 

123 transf = np.dot(transf, t) 

124 if self.debug >= 2: 

125 print("Rotating %f degrees around z axis." % (-theta * degree,)) 

126 print(self.basis) 

127 assert abs(np.linalg.det(transf) - 1) < 1e-6 

128 # Rotate second basis vector into xy plane 

129 theta = math.atan2(self.basis[1, 2], self.basis[1, 1]) 

130 t = np.array([[1, 0, 0], 

131 [0, np.cos(theta), -np.sin(theta)], 

132 [0, np.sin(theta), np.cos(theta)]]) 

133 self.basis = np.dot(self.basis, t) 

134 transf = np.dot(transf, t) 

135 if self.debug >= 2: 

136 print("Rotating %f degrees around x axis." % (-theta * degree,)) 

137 print(self.basis) 

138 assert abs(np.linalg.det(transf) - 1) < 1e-6 

139 # Now we better rotate the atoms as well 

140 self.atoms = np.dot(self.atoms, transf) 

141 # ... and rotate miller_basis 

142 self.miller_basis = np.dot(self.miller_basis, transf) 

143 

144 def make_list_of_atoms(self): 

145 "Repeat the unit cell." 

146 nrep = self.size[0] * self.size[1] * self.size[2] 

147 if nrep <= 0: 

148 raise ValueError( 

149 "Cannot create a non-positive number of unit cells") 

150 # Now the unit cells must be merged. 

151 a2 = [] 

152 e2 = [] 

153 for i in range(self.size[0]): 

154 offset = self.basis[0] * i 

155 a2.append(self.atoms + offset[np.newaxis, :]) 

156 e2.append(self.elements) 

157 atoms = np.concatenate(a2) 

158 elements = np.concatenate(e2) 

159 a2 = [] 

160 e2 = [] 

161 for j in range(self.size[1]): 

162 offset = self.basis[1] * j 

163 a2.append(atoms + offset[np.newaxis, :]) 

164 e2.append(elements) 

165 atoms = np.concatenate(a2) 

166 elements = np.concatenate(e2) 

167 a2 = [] 

168 e2 = [] 

169 for k in range(self.size[2]): 

170 offset = self.basis[2] * k 

171 a2.append(atoms + offset[np.newaxis, :]) 

172 e2.append(elements) 

173 atoms = np.concatenate(a2) 

174 elements = np.concatenate(e2) 

175 del a2, e2 

176 assert len(atoms) == nrep * len(self.atoms) 

177 basis = np.array([[self.size[0], 0, 0], 

178 [0, self.size[1], 0], 

179 [0, 0, self.size[2]]]) 

180 basis = np.dot(basis, self.basis) 

181 

182 # Tiny elements should be replaced by zero. The cutoff is 

183 # determined by chop_tolerance which is a class attribute. 

184 basis = np.where(np.abs(basis) < self.chop_tolerance, 

185 0.0, basis) 

186 

187 # None should be replaced, and memory should be freed. 

188 lattice = Lattice(positions=atoms, cell=basis, numbers=elements, 

189 pbc=self.pbc) 

190 lattice.millerbasis = self.miller_basis 

191 # Add info for lattice.surface.AddAdsorbate 

192 lattice._addsorbate_info_size = np.array(self.size[:2]) 

193 return lattice 

194 

195 def process_element(self, element): 

196 "Extract atomic number from element" 

197 # The types that can be elements: integers and strings 

198 if self.element_basis is None: 

199 if isinstance(element, str): 

200 self.atomicnumber = ase.data.atomic_numbers[element] 

201 elif isinstance(element, int): 

202 self.atomicnumber = element 

203 else: 

204 raise TypeError( 

205 "The symbol argument must be a string or an atomic number.") 

206 else: 

207 atomicnumber = [] 

208 try: 

209 if len(element) != max(self.element_basis) + 1: 

210 oops = True 

211 else: 

212 oops = False 

213 except TypeError: 

214 oops = True 

215 if oops: 

216 raise TypeError( 

217 ("The symbol argument must be a sequence of length %d" 

218 + " (one for each kind of lattice position") 

219 % (max(self.element_basis) + 1,)) 

220 for e in element: 

221 if isinstance(e, str): 

222 atomicnumber.append(ase.data.atomic_numbers[e]) 

223 elif isinstance(e, int): 

224 atomicnumber.append(e) 

225 else: 

226 raise TypeError( 

227 "The symbols argument must be a sequence of strings " 

228 "or atomic numbers.") 

229 self.atomicnumber = [atomicnumber[i] for i in self.element_basis] 

230 assert len(self.atomicnumber) == len(self.bravais_basis) 

231 

232 def convert_to_natural_basis(self): 

233 "Convert directions and miller indices to the natural basis." 

234 self.directions = np.dot(self.directions, self.inverse_basis) 

235 if self.bravais_basis is not None: 

236 self.natural_bravais_basis = np.dot(self.bravais_basis, 

237 self.inverse_basis) 

238 for i in (0, 1, 2): 

239 self.directions[i] = reduceindex(self.directions[i]) 

240 for i in (0, 1, 2): 

241 (j, k) = self.other[i] 

242 self.miller[i] = reduceindex(self.handedness * 

243 cross(self.directions[j], 

244 self.directions[k])) 

245 

246 def calc_num_atoms(self): 

247 v = int(round(abs(np.linalg.det(self.directions)))) 

248 if self.bravais_basis is None: 

249 return v 

250 else: 

251 return v * len(self.bravais_basis) 

252 

253 def make_unit_cell(self): 

254 "Make the unit cell." 

255 # Make three loops, and find the positions in the integral 

256 # lattice. Each time a position is found, the atom is placed 

257 # in the real unit cell by put_atom(). 

258 self.natoms = self.calc_num_atoms() 

259 self.nput = 0 

260 self.atoms = np.zeros((self.natoms, 3), float) 

261 self.elements = np.zeros(self.natoms, int) 

262 self.farpoint = sum(self.directions) 

263 # printprogress = self.debug and (len(self.atoms) > 250) 

264 # Find the radius of the sphere containing the whole system 

265 sqrad = 0 

266 for i in (0, 1): 

267 for j in (0, 1): 

268 for k in (0, 1): 

269 vect = (i * self.directions[0] + 

270 j * self.directions[1] + 

271 k * self.directions[2]) 

272 if np.dot(vect, vect) > sqrad: 

273 sqrad = np.dot(vect, vect) 

274 del i, j, k 

275 # Loop along first crystal axis (i) 

276 for (istart, istep) in ((0, 1), (-1, -1)): 

277 i = istart 

278 icont = True 

279 while icont: 

280 nj = 0 

281 for (jstart, jstep) in ((0, 1), (-1, -1)): 

282 j = jstart 

283 jcont = True 

284 while jcont: 

285 nk = 0 

286 for (kstart, kstep) in ((0, 1), (-1, -1)): 

287 k = kstart 

288 kcont = True 

289 while kcont: 

290 # Now (i,j,k) loops over Z^3, except that 

291 # the loops can be cut off when we get outside 

292 # the unit cell. 

293 point = np.array((i, j, k)) 

294 if self.inside(point): 

295 self.put_atom(point) 

296 nk += 1 

297 nj += 1 

298 # Is k too high? 

299 if np.dot(point, point) > sqrad: 

300 assert not self.inside(point) 

301 kcont = False 

302 k += kstep 

303 # Is j too high? 

304 if i * i + j * j > sqrad: 

305 jcont = False 

306 j += jstep 

307 # Is i too high? 

308 if i * i > sqrad: 

309 icont = False 

310 i += istep 

311 # if printprogress: 

312 # perce = int(100*self.nput / len(self.atoms)) 

313 # if perce > percent + 10: 

314 # print ("%d%%" % perce), 

315 # percent = perce 

316 assert(self.nput == self.natoms) 

317 

318 def inside(self, point): 

319 "Is a point inside the unit cell?" 

320 return (np.dot(self.miller[0], point) >= 0 and 

321 np.dot(self.miller[0], point - self.farpoint) < 0 and 

322 np.dot(self.miller[1], point) >= 0 and 

323 np.dot(self.miller[1], point - self.farpoint) < 0 and 

324 np.dot(self.miller[2], point) >= 0 and 

325 np.dot(self.miller[2], point - self.farpoint) < 0) 

326 

327 def put_atom(self, point): 

328 "Place an atom given its integer coordinates." 

329 if self.bravais_basis is None: 

330 # No basis - just place a single atom 

331 pos = np.dot(point, self.crystal_basis) 

332 if self.debug >= 2: 

333 print('Placing an atom at (%d,%d,%d) ~ (%.3f, %.3f, %.3f).' % 

334 (tuple(point) + tuple(pos))) 

335 self.atoms[self.nput] = pos 

336 self.elements[self.nput] = self.atomicnumber 

337 self.nput += 1 

338 else: 

339 for i, offset in enumerate(self.natural_bravais_basis): 

340 pos = np.dot(point + offset, self.crystal_basis) 

341 if self.debug >= 2: 

342 print('Placing an atom at (%d+%f, %d+%f, %d+%f) ~ ' 

343 '(%.3f, %.3f, %.3f).' % 

344 (point[0], offset[0], point[1], offset[1], 

345 point[2], offset[2], pos[0], pos[1], pos[2])) 

346 self.atoms[self.nput] = pos 

347 if self.element_basis is None: 

348 self.elements[self.nput] = self.atomicnumber 

349 else: 

350 self.elements[self.nput] = self.atomicnumber[i] 

351 self.nput += 1 

352 

353 def find_directions(self, directions, miller): 

354 """ 

355 Find missing directions and miller indices from the specified ones. 

356 """ 

357 directions = np.asarray(directions).tolist() 

358 miller = list(miller) # np.asarray(miller).tolist() 

359 # If no directions etc are specified, use a sensible default. 

360 if directions == [None, None, None] and miller == [None, None, None]: 

361 directions = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] 

362 # Now fill in missing directions and miller indices. This is an 

363 # iterative process. 

364 change = 1 

365 while change: 

366 change = False 

367 missing = 0 

368 for i in (0, 1, 2): 

369 j, k = self.other[i] 

370 if directions[i] is None: 

371 missing += 1 

372 if miller[j] is not None and miller[k] is not None: 

373 directions[i] = reduceindex(cross(miller[j], 

374 miller[k])) 

375 change = True 

376 if self.debug >= 2: 

377 print( 

378 "Calculating directions[%d] from miller " 

379 "indices" % i) 

380 if miller[i] is None: 

381 missing += 1 

382 if directions[j] is not None and directions[k] is not None: 

383 miller[i] = reduceindex(cross(directions[j], 

384 directions[k])) 

385 change = True 

386 if self.debug >= 2: 

387 print("Calculating miller[%d] from directions" % i) 

388 if missing: 

389 raise ValueError( 

390 "Specification of directions and miller indices is incomplete.") 

391 # Make sure that everything is Numeric arrays 

392 self.directions = np.array(directions) 

393 self.miller = np.array(miller) 

394 # Check for zero-volume unit cell 

395 if abs(np.linalg.det(self.directions)) < 1e-10: 

396 raise ValueError( 

397 "The direction vectors are linearly dependent " 

398 "(unit cell volume would be zero)") 

399 # Check for left-handed coordinate system 

400 if np.linalg.det(self.directions) < 0: 

401 print("WARNING: Creating a left-handed coordinate system!") 

402 self.miller = -self.miller 

403 self.handedness = -1 

404 else: 

405 self.handedness = 1 

406 # Now check for consistency 

407 for i in (0, 1, 2): 

408 (j, k) = self.other[i] 

409 m = reduceindex(self.handedness * 

410 cross(self.directions[j], self.directions[k])) 

411 if sum(np.not_equal(m, self.miller[i])): 

412 print( 

413 "ERROR: Miller index %s is inconsisten with " 

414 "directions %d and %d" % (i, j, k)) 

415 print("Miller indices:") 

416 print(str(self.miller)) 

417 print("Directions:") 

418 print(str(self.directions)) 

419 raise ValueError( 

420 "Inconsistent specification of miller indices and " 

421 "directions.") 

422 

423 def print_directions_and_miller(self, txt=""): 

424 "Print direction vectors and Miller indices." 

425 print("Direction vectors of unit cell%s:" % (txt,)) 

426 for i in (0, 1, 2): 

427 print(" ", self.directions[i]) 

428 print("Miller indices of surfaces%s:" % (txt,)) 

429 for i in (0, 1, 2): 

430 print(" ", self.miller[i]) 

431 

432 

433class MillerInfo: 

434 """Mixin class to provide information about Miller indices.""" 

435 

436 def miller_to_direction(self, miller): 

437 """Returns the direction corresponding to a given Miller index.""" 

438 return np.dot(miller, self.millerbasis) 

439 

440 

441class Lattice(Atoms, MillerInfo): 

442 """List of atoms initially containing a regular lattice of atoms. 

443 

444 A part from the usual list of atoms methods this list of atoms type 

445 also has a method, `miller_to_direction`, used to convert from Miller 

446 indices to directions in the coordinate system of the lattice. 

447 """ 

448 pass 

449 

450 

451# Helper functions 

452def cross(a, b): 

453 """The cross product of two vectors.""" 

454 return np.array((a[1] * b[2] - b[1] * a[2], 

455 a[2] * b[0] - b[2] * a[0], 

456 a[0] * b[1] - b[0] * a[1])) 

457 

458 

459def reduceindex(M): 

460 """Reduce Miller index to the lowest equivalent integers.""" 

461 oldM = M 

462 g = math.gcd(M[0], M[1]) 

463 h = math.gcd(g, M[2]) 

464 while h != 1: 

465 if h == 0: 

466 raise ValueError( 

467 "Division by zero: Are the miller indices linearly dependent?") 

468 M = M // h 

469 g = math.gcd(M[0], M[1]) 

470 h = math.gcd(g, M[2]) 

471 if np.dot(oldM, M) > 0: 

472 return M 

473 else: 

474 return -M