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

1# fmt: off

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

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

6"""

8import math

9from typing import Optional, Sequence

11import numpy as np

13import ase.data

14from ase.atoms import Atoms

17class Bravais:

18 """Bravais lattice factory.

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

21 (SC, FCC, ...).

22 """

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

25 # in classes inhering from Bravais:

26 # get_lattice_constant

27 # make_crystal_basis

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

29 # in classes inhering from Bravais:

30 # int_basis

31 # inverse_basis

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

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

36 # None if no basis is present.

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

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

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

41 # bravais_basis could contain four vectors, and element_basis

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

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

44 # same type.

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

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

48 chop_tolerance = 1e-10

50 def __call__(self, symbol,

51 directions=(None, None, None), miller=(None, None, None),

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

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

54 "Create a lattice."

55 self.size = size

56 self.pbc = pbc

57 self.debug = debug

58 self.process_element(symbol)

59 self.find_directions(directions, miller)

60 if self.debug:

61 self.print_directions_and_miller()

62 self.convert_to_natural_basis()

63 if self.debug >= 2:

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

65 if latticeconstant is None:

66 if self.element_basis is None:

67 self.latticeconstant = self.get_lattice_constant()

68 else:

69 raise ValueError(

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

71 else:

72 self.latticeconstant = latticeconstant

73 if self.debug:

74 print(

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

76 self.calc_num_atoms())

77 if self.debug >= 2:

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

79 if self.bravais_basis is not None:

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

81 self.make_crystal_basis()

82 self.make_unit_cell()

83 if align:

84 self.align()

85 return self.make_list_of_atoms()

87 def align(self):

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

89 degree = 180 / np.pi

90 if self.debug >= 2:

91 print("Basis before alignment:")

92 print(self.basis)

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

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

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

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

97 [1, 0, 0],

98 [0, 0, 1]], float)

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

100 transf = t

101 if self.debug >= 2:

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

103 "stability.")

104 print(self.basis)

105 else:

106 transf = np.identity(3, float)

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

108 # Rotate first basis vector into xy plane

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

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

111 [0, 1, 0],

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

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

114 transf = np.dot(transf, t)

115 if self.debug >= 2:

116 print(f"Rotating {-theta * degree:f} degrees around y axis.")

117 print(self.basis)

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

119 # Rotate first basis vector to point along x axis

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

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

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

123 [0, 0, 1]])

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

125 transf = np.dot(transf, t)

126 if self.debug >= 2:

127 print(f"Rotating {-theta * degree:f} degrees around z axis.")

128 print(self.basis)

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

130 # Rotate second basis vector into xy plane

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

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

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

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

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

136 transf = np.dot(transf, t)

137 if self.debug >= 2:

138 print(f"Rotating {-theta * degree:f} degrees around x axis.")

139 print(self.basis)

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

141 # Now we better rotate the atoms as well

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

143 # ... and rotate miller_basis

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

145

146 def make_list_of_atoms(self):

147 "Repeat the unit cell."

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

149 if nrep <= 0:

150 raise ValueError(

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

152 # Now the unit cells must be merged.

153 a2 = []

154 e2 = []

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

156 offset = self.basis[0] * i

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

158 e2.append(self.elements)

159 atoms = np.concatenate(a2)

160 elements = np.concatenate(e2)

161 a2 = []

162 e2 = []

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

164 offset = self.basis[1] * j

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

166 e2.append(elements)

167 atoms = np.concatenate(a2)

168 elements = np.concatenate(e2)

169 a2 = []

170 e2 = []

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

172 offset = self.basis[2] * k

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

174 e2.append(elements)

175 atoms = np.concatenate(a2)

176 elements = np.concatenate(e2)

177 del a2, e2

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

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

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

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

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

183

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

185 # determined by chop_tolerance which is a class attribute.

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

187 0.0, basis)

188

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

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

191 pbc=self.pbc)

192 lattice.millerbasis = self.miller_basis

193 # Add info for lattice.surface.AddAdsorbate

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

195 return lattice

196

197 def process_element(self, element):

198 "Extract atomic number from element"

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

200 if self.element_basis is None:

201 if isinstance(element, str):

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

203 elif isinstance(element, int):

204 self.atomicnumber = element

205 else:

206 raise TypeError(

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

208 else:

209 atomicnumber = []

210 try:

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

212 oops = True

213 else:

214 oops = False

215 except TypeError:

216 oops = True

217 if oops:

218 raise TypeError(

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

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

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

222 for e in element:

223 if isinstance(e, str):

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

225 elif isinstance(e, int):

226 atomicnumber.append(e)

227 else:

228 raise TypeError(

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

230 "or atomic numbers.")

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

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

233

234 def convert_to_natural_basis(self):

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

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

237 if self.bravais_basis is not None:

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

239 self.inverse_basis)

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

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

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

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

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

245 cross(self.directions[j],

246 self.directions[k]))

247

248 def calc_num_atoms(self):

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

250 if self.bravais_basis is None:

251 return v

252 else:

253 return v * len(self.bravais_basis)

254

255 def make_unit_cell(self):

256 "Make the unit cell."

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

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

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

260 self.natoms = self.calc_num_atoms()

261 self.nput = 0

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

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

264 self.farpoint = sum(self.directions)

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

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

267 sqrad = 0

268 for i in (0, 1):

269 for j in (0, 1):

270 for k in (0, 1):

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

272 j * self.directions[1] +

273 k * self.directions[2])

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

275 sqrad = np.dot(vect, vect)

276 del i, j, k

277 # Loop along first crystal axis (i)

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

279 i = istart

280 icont = True

281 while icont:

282 nj = 0

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

284 j = jstart

285 jcont = True

286 while jcont:

287 nk = 0

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

289 k = kstart

290 kcont = True

291 while kcont:

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

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

294 # the unit cell.

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

296 if self.inside(point):

297 self.put_atom(point)

298 nk += 1

299 nj += 1

300 # Is k too high?

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

302 assert not self.inside(point)

303 kcont = False

304 k += kstep

305 # Is j too high?

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

307 jcont = False

308 j += jstep

309 # Is i too high?

310 if i * i > sqrad:

311 icont = False

312 i += istep

313 # if printprogress:

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

315 # if perce > percent + 10:

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

317 # percent = perce

318 assert self.nput == self.natoms

319

320 def inside(self, point):

321 "Is a point inside the unit cell?"

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

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

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

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

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

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

328

329 def put_atom(self, point):

330 "Place an atom given its integer coordinates."

331 if self.bravais_basis is None:

332 # No basis - just place a single atom

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

334 if self.debug >= 2:

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

336 (tuple(point) + tuple(pos)))

337 self.atoms[self.nput] = pos

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

339 self.nput += 1

340 else:

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

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

343 if self.debug >= 2:

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

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

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

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

348 self.atoms[self.nput] = pos

349 if self.element_basis is None:

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

351 else:

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

353 self.nput += 1

354

355 def find_directions(self, directions, miller):

356 """

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

358 """

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

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

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

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

363 directions = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

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

365 # iterative process.

366 change = 1

367 while change:

368 change = False

369 missing = 0

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

371 j, k = self.other[i]

372 if directions[i] is None:

373 missing += 1

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

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

376 miller[k]))

377 change = True

378 if self.debug >= 2:

379 print(

380 "Calculating directions[%d] from miller "

381 "indices" % i)

382 if miller[i] is None:

383 missing += 1

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

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

386 directions[k]))

387 change = True

388 if self.debug >= 2:

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

390 if missing:

391 raise ValueError(

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

393 # Make sure that everything is Numeric arrays

394 self.directions = np.array(directions)

395 self.miller = np.array(miller)

396 # Check for zero-volume unit cell

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

398 raise ValueError(

399 "The direction vectors are linearly dependent "

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

401 # Check for left-handed coordinate system

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

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

404 self.miller = -self.miller

405 self.handedness = -1

406 else:

407 self.handedness = 1

408 # Now check for consistency

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

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

411 m = reduceindex(self.handedness *

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

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

414 print(

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

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

417 print("Miller indices:")

418 print(str(self.miller))

419 print("Directions:")

420 print(str(self.directions))

421 raise ValueError(

422 "Inconsistent specification of miller indices and "

423 "directions.")

424

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

426 "Print direction vectors and Miller indices."

427 print(f"Direction vectors of unit cell{txt}:")

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

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

430 print(f"Miller indices of surfaces{txt}:")

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

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

433

434

435class MillerInfo:

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

437

438 def miller_to_direction(self, miller):

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

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

441

442

443class Lattice(Atoms, MillerInfo):

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

445

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

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

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

449 """

450

451

452# Helper functions

453def cross(a, b):

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

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

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

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

458

459

460def reduceindex(M):

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

462 oldM = M

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

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

465 while h != 1:

466 if h == 0:

467 raise ValueError(

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

469 M = M // h

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

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

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

473 return M

474 else:

475 return -M