Coverage for /builds/kinetik161/ase/ase/lattice/bravais.py: 81.18%

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

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

4"""

6import math

7from typing import Optional, Sequence

9import numpy as np

11import ase.data

12from ase.atoms import Atoms

15class Bravais:

16 """Bravais lattice factory.

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

19 (SC, FCC, ...).

20 """

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

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

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

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

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

46 chop_tolerance = 1e-10

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()

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(f"Rotating {-theta * degree:f} degrees around y axis.")

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(f"Rotating {-theta * degree:f} degrees around z axis.")

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(f"Rotating {-theta * degree:f} degrees around x axis.")

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(f"Direction vectors of unit cell{txt}:")

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

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

428 print(f"Miller indices of surfaces{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

449

450# Helper functions

451def cross(a, b):

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

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

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

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

456

457

458def reduceindex(M):

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

460 oldM = M

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

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

463 while h != 1:

464 if h == 0:

465 raise ValueError(

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

467 M = M // h

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

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

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

471 return M

472 else:

473 return -M