Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import re 

2import warnings 

3from typing import Dict 

4 

5import numpy as np 

6 

7import ase # Annotations 

8from ase.utils import jsonable 

9from ase.cell import Cell 

10 

11 

12def monkhorst_pack(size): 

13 """Construct a uniform sampling of k-space of given size.""" 

14 if np.less_equal(size, 0).any(): 

15 raise ValueError('Illegal size: %s' % list(size)) 

16 kpts = np.indices(size).transpose((1, 2, 3, 0)).reshape((-1, 3)) 

17 return (kpts + 0.5) / size - 0.5 

18 

19 

20def get_monkhorst_pack_size_and_offset(kpts): 

21 """Find Monkhorst-Pack size and offset. 

22 

23 Returns (size, offset), where:: 

24 

25 kpts = monkhorst_pack(size) + offset. 

26 

27 The set of k-points must not have been symmetry reduced.""" 

28 

29 if len(kpts) == 1: 

30 return np.ones(3, int), np.array(kpts[0], dtype=float) 

31 

32 size = np.zeros(3, int) 

33 for c in range(3): 

34 # Determine increment between k-points along current axis 

35 delta = max(np.diff(np.sort(kpts[:, c]))) 

36 

37 # Determine number of k-points as inverse of distance between kpoints 

38 if delta > 1e-8: 

39 size[c] = int(round(1.0 / delta)) 

40 else: 

41 size[c] = 1 

42 

43 if size.prod() == len(kpts): 

44 kpts0 = monkhorst_pack(size) 

45 offsets = kpts - kpts0 

46 

47 # All offsets must be identical: 

48 if (offsets.ptp(axis=0) < 1e-9).all(): 

49 return size, offsets[0].copy() 

50 

51 raise ValueError('Not an ASE-style Monkhorst-Pack grid!') 

52 

53 

54def mindistance2monkhorstpack(atoms, *, 

55 min_distance, 

56 maxperdim=16, 

57 even=True): 

58 """Find a Monkhorst-Pack grid (nx, ny, nz) with lowest number of 

59 k-points in the *reducible* Brillouin zone, which still satisfying 

60 a given minimum distance (`min_distance`) condition in real space 

61 (nx, ny, nz)-supercell. 

62 

63 Compared to ase.calculators.calculator kptdensity2monkhorstpack 

64 routine, this metric is based on a physical quantity (real space 

65 distance), and it doesn't depend on non-physical quantities, such as 

66 the cell vectors, since basis vectors can be always transformed 

67 with integer determinant one matrices. In other words, it is 

68 invariant to particular choice of cell representations. 

69 

70 On orthogonal cells, min_distance = 2 * np.pi * kptdensity. 

71 """ 

72 return _mindistance2monkhorstpack(atoms.cell, atoms.pbc, 

73 min_distance, maxperdim, even) 

74 

75 

76def _mindistance2monkhorstpack(cell, pbc_c, min_distance, maxperdim, even): 

77 from ase import Atoms 

78 from ase.neighborlist import NeighborList 

79 

80 step = 2 if even else 1 

81 nl = NeighborList([min_distance / 2], skin=0.0, 

82 self_interaction=False, bothways=False) 

83 

84 def check(nkpts_c): 

85 nl.update(Atoms('H', cell=cell @ np.diag(nkpts_c), pbc=pbc_c)) 

86 return len(nl.get_neighbors(0)[1]) == 0 

87 

88 def generate_mpgrids(): 

89 ranges = [range(step, maxperdim + 1, step) 

90 if pbc else range(1, 2) for pbc in pbc_c] 

91 nkpts_nc = np.column_stack([*map(np.ravel, np.meshgrid(*ranges))]) 

92 yield from sorted(nkpts_nc, key=lambda nkpts_c: np.prod(nkpts_c)) 

93 

94 try: 

95 return next(filter(check, generate_mpgrids())) 

96 except StopIteration: 

97 raise ValueError('Could not find a proper k-point grid for the system.' 

98 ' Try running with a larger maxperdim.') 

99 

100 

101def get_monkhorst_shape(kpts): 

102 warnings.warn('Use get_monkhorst_pack_size_and_offset()[0] instead.') 

103 return get_monkhorst_pack_size_and_offset(kpts)[0] 

104 

105 

106def kpoint_convert(cell_cv, skpts_kc=None, ckpts_kv=None): 

107 """Convert k-points between scaled and cartesian coordinates. 

108 

109 Given the atomic unit cell, and either the scaled or cartesian k-point 

110 coordinates, the other is determined. 

111 

112 The k-point arrays can be either a single point, or a list of points, 

113 i.e. the dimension k can be empty or multidimensional. 

114 """ 

115 if ckpts_kv is None: 

116 icell_cv = 2 * np.pi * np.linalg.pinv(cell_cv).T 

117 return np.dot(skpts_kc, icell_cv) 

118 elif skpts_kc is None: 

119 return np.dot(ckpts_kv, cell_cv.T) / (2 * np.pi) 

120 else: 

121 raise KeyError('Either scaled or cartesian coordinates must be given.') 

122 

123 

124def parse_path_string(s): 

125 """Parse compact string representation of BZ path. 

126 

127 A path string can have several non-connected sections separated by 

128 commas. The return value is a list of sections where each section is a 

129 list of labels. 

130 

131 Examples: 

132 

133 >>> parse_path_string('GX') 

134 [['G', 'X']] 

135 >>> parse_path_string('GX,M1A') 

136 [['G', 'X'], ['M1', 'A']] 

137 """ 

138 paths = [] 

139 for path in s.split(','): 

140 names = [name 

141 for name in re.split(r'([A-Z][a-z0-9]*)', path) 

142 if name] 

143 paths.append(names) 

144 return paths 

145 

146 

147def resolve_kpt_path_string(path, special_points): 

148 paths = parse_path_string(path) 

149 coords = [np.array([special_points[sym] for sym in subpath]).reshape(-1, 3) 

150 for subpath in paths] 

151 return paths, coords 

152 

153 

154def resolve_custom_points(pathspec, special_points, eps): 

155 """Resolve a path specification into a string. 

156 

157 The path specification is a list path segments, each segment being a kpoint 

158 label or kpoint coordinate, or a single such segment. 

159 

160 Return a string representing the same path. Generic kpoint labels 

161 are generated dynamically as necessary, updating the special_point 

162 dictionary if necessary. The tolerance eps is used to see whether 

163 coordinates are close enough to a special point to deserve being 

164 labelled as such.""" 

165 # This should really run on Cartesian coordinates but we'll probably 

166 # be lazy and call it on scaled ones. 

167 

168 # We may add new points below so take a copy of the input: 

169 special_points = dict(special_points) 

170 

171 if len(pathspec) == 0: 

172 return '', special_points 

173 

174 if isinstance(pathspec, str): 

175 pathspec = parse_path_string(pathspec) 

176 

177 def looks_like_single_kpoint(obj): 

178 if isinstance(obj, str): 

179 return True 

180 try: 

181 arr = np.asarray(obj, float) 

182 except ValueError: 

183 return False 

184 else: 

185 return arr.shape == (3,) 

186 

187 # We accept inputs that are either 

188 # 1) string notation 

189 # 2) list of kpoints (each either a string or three floats) 

190 # 3) list of list of kpoints; each toplevel list is a contiguous subpath 

191 # Here we detect form 2 and normalize to form 3: 

192 for thing in pathspec: 

193 if looks_like_single_kpoint(thing): 

194 pathspec = [pathspec] 

195 break 

196 

197 def name_generator(): 

198 counter = 0 

199 while True: 

200 name = 'Kpt{}'.format(counter) 

201 yield name 

202 counter += 1 

203 custom_names = name_generator() 

204 

205 labelseq = [] 

206 for subpath in pathspec: 

207 for kpt in subpath: 

208 if isinstance(kpt, str): 

209 if kpt not in special_points: 

210 raise KeyError('No kpoint "{}" among "{}"' 

211 .format(kpt, 

212 ''.join(special_points))) 

213 labelseq.append(kpt) 

214 continue 

215 

216 kpt = np.asarray(kpt, float) 

217 if not kpt.shape == (3,): 

218 raise ValueError(f'Not a valid kpoint: {kpt}') 

219 

220 for key, val in special_points.items(): 

221 if np.abs(kpt - val).max() < eps: 

222 labelseq.append(key) 

223 break # Already present 

224 else: 

225 # New special point - search for name we haven't used yet: 

226 name = next(custom_names) 

227 while name in special_points: 

228 name = next(custom_names) 

229 special_points[name] = kpt 

230 labelseq.append(name) 

231 labelseq.append(',') 

232 

233 last = labelseq.pop() 

234 assert last == ',' 

235 return ''.join(labelseq), special_points 

236 

237 

238def normalize_special_points(special_points): 

239 dct = {} 

240 for name, value in special_points.items(): 

241 if not isinstance(name, str): 

242 raise TypeError('Expected name to be a string') 

243 if not np.shape(value) == (3,): 

244 raise ValueError('Expected 3 kpoint coordinates') 

245 dct[name] = np.asarray(value, float) 

246 return dct 

247 

248 

249@jsonable('bandpath') 

250class BandPath: 

251 """Represents a Brillouin zone path or bandpath. 

252 

253 A band path has a unit cell, a path specification, special points, 

254 and interpolated k-points. Band paths are typically created 

255 indirectly using the :class:`~ase.geometry.Cell` or 

256 :class:`~ase.lattice.BravaisLattice` classes: 

257 

258 >>> from ase.lattice import CUB 

259 >>> path = CUB(3).bandpath() 

260 >>> path 

261 BandPath(path='GXMGRX,MR', cell=[3x3], special_points={GMRX}, kpts=[40x3]) 

262 

263 Band paths support JSON I/O: 

264 

265 >>> from ase.io.jsonio import read_json 

266 >>> path.write('mybandpath.json') 

267 >>> read_json('mybandpath.json') 

268 BandPath(path='GXMGRX,MR', cell=[3x3], special_points={GMRX}, kpts=[40x3]) 

269 

270 """ 

271 def __init__(self, cell, kpts=None, 

272 special_points=None, path=None): 

273 if kpts is None: 

274 kpts = np.empty((0, 3)) 

275 

276 if special_points is None: 

277 special_points = {} 

278 else: 

279 special_points = normalize_special_points(special_points) 

280 

281 if path is None: 

282 path = '' 

283 

284 cell = Cell(cell) 

285 self._cell = cell 

286 kpts = np.asarray(kpts) 

287 assert kpts.ndim == 2 and kpts.shape[1] == 3 and kpts.dtype == float 

288 self._icell = self.cell.reciprocal() 

289 self._kpts = kpts 

290 self._special_points = special_points 

291 if not isinstance(path, str): 

292 raise TypeError(f'path must be a string; was {path!r}') 

293 self._path = path 

294 

295 @property 

296 def cell(self) -> Cell: 

297 """The :class:`~ase.cell.Cell` of this BandPath.""" 

298 return self._cell 

299 

300 @property 

301 def icell(self) -> Cell: 

302 """Reciprocal cell of this BandPath as a :class:`~ase.cell.Cell`.""" 

303 return self._icell 

304 

305 @property 

306 def kpts(self) -> np.ndarray: 

307 """The kpoints of this BandPath as an array of shape (nkpts, 3). 

308 

309 The kpoints are given in units of the reciprocal cell.""" 

310 return self._kpts 

311 

312 @property 

313 def special_points(self) -> Dict[str, np.ndarray]: 

314 """Special points of this BandPath as a dictionary. 

315 

316 The dictionary maps names (such as `'G'`) to kpoint coordinates 

317 in units of the reciprocal cell as a 3-element numpy array. 

318 

319 It's unwise to edit this dictionary directly. If you need that, 

320 consider deepcopying it.""" 

321 return self._special_points 

322 

323 @property 

324 def path(self) -> str: 

325 """The string specification of this band path. 

326 

327 This is a specification of the form `'GXWKGLUWLK,UX'`. 

328 

329 Comma marks a discontinuous jump: K is not connected to U.""" 

330 return self._path 

331 

332 def transform(self, op: np.ndarray) -> 'BandPath': 

333 """Apply 3x3 matrix to this BandPath and return new BandPath. 

334 

335 This is useful for converting the band path to another cell. 

336 The operation will typically be a permutation/flipping 

337 established by a function such as Niggli reduction.""" 

338 # XXX acceptable operations are probably only those 

339 # who come from Niggli reductions (permutations etc.) 

340 # 

341 # We should insert a check. 

342 # I wonder which operations are valid? They won't be valid 

343 # if they change lengths, volume etc. 

344 special_points = {} 

345 for name, value in self.special_points.items(): 

346 special_points[name] = value @ op 

347 

348 return BandPath(op.T @ self.cell, kpts=self.kpts @ op, 

349 special_points=special_points, 

350 path=self.path) 

351 

352 def todict(self): 

353 return {'kpts': self.kpts, 

354 'special_points': self.special_points, 

355 'labelseq': self.path, 

356 'cell': self.cell} 

357 

358 def interpolate( 

359 self, 

360 path: str = None, 

361 npoints: int = None, 

362 special_points: Dict[str, np.ndarray] = None, 

363 density: float = None, 

364 ) -> 'BandPath': 

365 """Create new bandpath, (re-)interpolating kpoints from this one.""" 

366 if path is None: 

367 path = self.path 

368 

369 if special_points is None: 

370 special_points = self.special_points 

371 

372 pathnames, pathcoords = resolve_kpt_path_string(path, special_points) 

373 kpts, x, X = paths2kpts(pathcoords, self.cell, npoints, density) 

374 return BandPath(self.cell, kpts, path=path, 

375 special_points=special_points) 

376 

377 def _scale(self, coords): 

378 return np.dot(coords, self.icell) 

379 

380 def __repr__(self): 

381 return ('{}(path={}, cell=[3x3], special_points={{{}}}, kpts=[{}x3])' 

382 .format(self.__class__.__name__, 

383 repr(self.path), 

384 ''.join(sorted(self.special_points)), 

385 len(self.kpts))) 

386 

387 def cartesian_kpts(self) -> np.ndarray: 

388 """Get Cartesian kpoints from this bandpath.""" 

389 return self._scale(self.kpts) 

390 

391 def __iter__(self): 

392 """XXX Compatibility hack for bandpath() function. 

393 

394 bandpath() now returns a BandPath object, which is a Good 

395 Thing. However it used to return a tuple of (kpts, x_axis, 

396 special_x_coords), and people would use tuple unpacking for 

397 those. 

398 

399 This function makes tuple unpacking work in the same way. 

400 It will be removed in the future. 

401 

402 """ 

403 import warnings 

404 warnings.warn('Please do not use (kpts, x, X) = bandpath(...). ' 

405 'Use path = bandpath(...) and then kpts = path.kpts and ' 

406 '(x, X, labels) = path.get_linear_kpoint_axis().') 

407 yield self.kpts 

408 

409 x, xspecial, _ = self.get_linear_kpoint_axis() 

410 yield x 

411 yield xspecial 

412 

413 def __getitem__(self, index): 

414 # Temp compatibility stuff, see __iter__ 

415 return tuple(self)[index] 

416 

417 def get_linear_kpoint_axis(self, eps=1e-5): 

418 """Define x axis suitable for plotting a band structure. 

419 

420 See :func:`ase.dft.kpoints.labels_from_kpts`.""" 

421 

422 index2name = self._find_special_point_indices(eps) 

423 indices = sorted(index2name) 

424 labels = [index2name[index] for index in indices] 

425 xcoords, special_xcoords = indices_to_axis_coords( 

426 indices, self.kpts, self.cell) 

427 return xcoords, special_xcoords, labels 

428 

429 def _find_special_point_indices(self, eps): 

430 """Find indices of kpoints which are close to special points. 

431 

432 The result is returned as a dictionary mapping indices to labels.""" 

433 # XXX must work in Cartesian coordinates for comparison to eps 

434 # to fully make sense! 

435 index2name = {} 

436 nkpts = len(self.kpts) 

437 

438 for name, kpt in self.special_points.items(): 

439 displacements = self.kpts - kpt[np.newaxis, :] 

440 distances = np.linalg.norm(displacements, axis=1) 

441 args = np.argwhere(distances < eps) 

442 for arg in args.flat: 

443 dist = distances[arg] 

444 # Check if an adjacent point exists and is even closer: 

445 neighbours = distances[max(arg - 1, 0):min(arg + 1, nkpts - 1)] 

446 if not any(neighbours < dist): 

447 index2name[arg] = name 

448 

449 return index2name 

450 

451 def plot(self, **plotkwargs): 

452 """Visualize this bandpath. 

453 

454 Plots the irreducible Brillouin zone and this bandpath.""" 

455 import ase.dft.bz as bz 

456 

457 # We previously had a "dimension=3" argument which is now unused. 

458 plotkwargs.pop('dimension', None) 

459 

460 special_points = self.special_points 

461 labelseq, coords = resolve_kpt_path_string(self.path, 

462 special_points) 

463 

464 paths = [] 

465 points_already_plotted = set() 

466 for subpath_labels, subpath_coords in zip(labelseq, coords): 

467 subpath_coords = np.array(subpath_coords) 

468 points_already_plotted.update(subpath_labels) 

469 paths.append((subpath_labels, self._scale(subpath_coords))) 

470 

471 # Add each special point as a single-point subpath if they were 

472 # not plotted already: 

473 for label, point in special_points.items(): 

474 if label not in points_already_plotted: 

475 paths.append(([label], [self._scale(point)])) 

476 

477 kw = {'vectors': True, 

478 'pointstyle': {'marker': '.'}} 

479 

480 kw.update(plotkwargs) 

481 return bz.bz_plot(self.cell, paths=paths, 

482 points=self.cartesian_kpts(), 

483 **kw) 

484 

485 def free_electron_band_structure( 

486 self, **kwargs 

487 ) -> 'ase.spectrum.band_structure.BandStructure': 

488 """Return band structure of free electrons for this bandpath. 

489 

490 Keyword arguments are passed to 

491 :class:`~ase.calculators.test.FreeElectrons`. 

492 

493 This is for mostly testing and visualization.""" 

494 from ase import Atoms 

495 from ase.calculators.test import FreeElectrons 

496 from ase.spectrum.band_structure import calculate_band_structure 

497 atoms = Atoms(cell=self.cell, pbc=True) 

498 atoms.calc = FreeElectrons(**kwargs) 

499 bs = calculate_band_structure(atoms, path=self) 

500 return bs 

501 

502 

503def bandpath(path, cell, npoints=None, density=None, special_points=None, 

504 eps=2e-4): 

505 """Make a list of kpoints defining the path between the given points. 

506 

507 path: list or str 

508 Can be: 

509 

510 * a string that parse_path_string() understands: 'GXL' 

511 * a list of BZ points: [(0, 0, 0), (0.5, 0, 0)] 

512 * or several lists of BZ points if the the path is not continuous. 

513 cell: 3x3 

514 Unit cell of the atoms. 

515 npoints: int 

516 Length of the output kpts list. If too small, at least the beginning 

517 and ending point of each path segment will be used. If None (default), 

518 it will be calculated using the supplied density or a default one. 

519 density: float 

520 k-points per 1/A on the output kpts list. If npoints is None, 

521 the number of k-points in the output list will be: 

522 npoints = density * path total length (in Angstroms). 

523 If density is None (default), use 5 k-points per A⁻¹. 

524 If the calculated npoints value is less than 50, a minimum value of 50 

525 will be used. 

526 special_points: dict or None 

527 Dictionary mapping names to special points. If None, the special 

528 points will be derived from the cell. 

529 eps: float 

530 Precision used to identify Bravais lattice, deducing special points. 

531 

532 You may define npoints or density but not both. 

533 

534 Return a :class:`~ase.dft.kpoints.BandPath` object.""" 

535 

536 cell = Cell.ascell(cell) 

537 return cell.bandpath(path, npoints=npoints, density=density, 

538 special_points=special_points, eps=eps) 

539 

540 

541DEFAULT_KPTS_DENSITY = 5 # points per 1/Angstrom 

542 

543 

544def paths2kpts(paths, cell, npoints=None, density=None): 

545 if not (npoints is None or density is None): 

546 raise ValueError('You may define npoints or density, but not both.') 

547 points = np.concatenate(paths) 

548 dists = points[1:] - points[:-1] 

549 lengths = [np.linalg.norm(d) for d in kpoint_convert(cell, skpts_kc=dists)] 

550 

551 i = 0 

552 for path in paths[:-1]: 

553 i += len(path) 

554 lengths[i - 1] = 0 

555 

556 length = sum(lengths) 

557 

558 if npoints is None: 

559 if density is None: 

560 density = DEFAULT_KPTS_DENSITY 

561 # Set npoints using the length of the path 

562 npoints = int(round(length * density)) 

563 

564 kpts = [] 

565 x0 = 0 

566 x = [] 

567 X = [0] 

568 for P, d, L in zip(points[:-1], dists, lengths): 

569 diff = length - x0 

570 if abs(diff) < 1e-6: 

571 n = 0 

572 else: 

573 n = max(2, int(round(L * (npoints - len(x)) / diff))) 

574 

575 for t in np.linspace(0, 1, n)[:-1]: 

576 kpts.append(P + t * d) 

577 x.append(x0 + t * L) 

578 x0 += L 

579 X.append(x0) 

580 if len(points): 

581 kpts.append(points[-1]) 

582 x.append(x0) 

583 

584 if len(kpts) == 0: 

585 kpts = np.empty((0, 3)) 

586 

587 return np.array(kpts), np.array(x), np.array(X) 

588 

589 

590get_bandpath = bandpath # old name 

591 

592 

593def find_bandpath_kinks(cell, kpts, eps=1e-5): 

594 """Find indices of those kpoints that are not interior to a line segment.""" 

595 # XXX Should use the Cartesian kpoints. 

596 # Else comparison to eps will be anisotropic and hence arbitrary. 

597 diffs = kpts[1:] - kpts[:-1] 

598 kinks = abs(diffs[1:] - diffs[:-1]).sum(1) > eps 

599 N = len(kpts) 

600 indices = [] 

601 if N > 0: 

602 indices.append(0) 

603 indices.extend(np.arange(1, N - 1)[kinks]) 

604 indices.append(N - 1) 

605 return indices 

606 

607 

608def labels_from_kpts(kpts, cell, eps=1e-5, special_points=None): 

609 """Get an x-axis to be used when plotting a band structure. 

610 

611 The first of the returned lists can be used as a x-axis when plotting 

612 the band structure. The second list can be used as xticks, and the third 

613 as xticklabels. 

614 

615 Parameters: 

616 

617 kpts: list 

618 List of scaled k-points. 

619 

620 cell: list 

621 Unit cell of the atomic structure. 

622 

623 Returns: 

624 

625 Three arrays; the first is a list of cumulative distances between k-points, 

626 the second is x coordinates of the special points, 

627 the third is the special points as strings. 

628 """ 

629 

630 if special_points is None: 

631 special_points = get_special_points(cell) 

632 points = np.asarray(kpts) 

633 # XXX Due to this mechanism, we are blind to special points 

634 # that lie on straight segments such as [K, G, -K]. 

635 indices = find_bandpath_kinks(cell, kpts, eps=1e-5) 

636 

637 labels = [] 

638 for kpt in points[indices]: 

639 for label, k in special_points.items(): 

640 if abs(kpt - k).sum() < eps: 

641 break 

642 else: 

643 # No exact match. Try modulus 1: 

644 for label, k in special_points.items(): 

645 if abs((kpt - k) % 1).sum() < eps: 

646 break 

647 else: 

648 label = '?' 

649 labels.append(label) 

650 

651 xcoords, ixcoords = indices_to_axis_coords(indices, points, cell) 

652 return xcoords, ixcoords, labels 

653 

654 

655def indices_to_axis_coords(indices, points, cell): 

656 jump = False # marks a discontinuity in the path 

657 xcoords = [0] 

658 for i1, i2 in zip(indices[:-1], indices[1:]): 

659 if not jump and i1 + 1 == i2: 

660 length = 0 

661 jump = True # we don't want two jumps in a row 

662 else: 

663 diff = points[i2] - points[i1] 

664 length = np.linalg.norm(kpoint_convert(cell, skpts_kc=diff)) 

665 jump = False 

666 xcoords.extend(np.linspace(0, length, i2 - i1 + 1)[1:] + xcoords[-1]) 

667 

668 xcoords = np.array(xcoords) 

669 return xcoords, xcoords[indices] 

670 

671 

672special_paths = { 

673 'cubic': 'GXMGRX,MR', 

674 'fcc': 'GXWKGLUWLK,UX', 

675 'bcc': 'GHNGPH,PN', 

676 'tetragonal': 'GXMGZRAZXR,MA', 

677 'orthorhombic': 'GXSYGZURTZ,YT,UX,SR', 

678 'hexagonal': 'GMKGALHA,LM,KH', 

679 'monoclinic': 'GYHCEM1AXH1,MDZ,YD', 

680 'rhombohedral type 1': 'GLB1,BZGX,QFP1Z,LP', 

681 'rhombohedral type 2': 'GPZQGFP1Q1LZ'} 

682 

683 

684def get_special_points(cell, lattice=None, eps=2e-4): 

685 """Return dict of special points. 

686 

687 The definitions are from a paper by Wahyu Setyawana and Stefano 

688 Curtarolo:: 

689 

690 https://doi.org/10.1016/j.commatsci.2010.05.010 

691 

692 cell: 3x3 ndarray 

693 Unit cell. 

694 lattice: str 

695 Optionally check that the cell is one of the following: cubic, fcc, 

696 bcc, orthorhombic, tetragonal, hexagonal or monoclinic. 

697 eps: float 

698 Tolerance for cell-check. 

699 """ 

700 

701 if isinstance(cell, str): 

702 warnings.warn('Please call this function with cell as the first ' 

703 'argument') 

704 lattice, cell = cell, lattice 

705 

706 cell = Cell.ascell(cell) 

707 # We create the bandpath because we want to transform the kpoints too, 

708 # from the canonical cell to the given one. 

709 # 

710 # Note that this function is missing a tolerance, epsilon. 

711 path = cell.bandpath(npoints=0) 

712 return path.special_points 

713 

714 

715def monkhorst_pack_interpolate(path, values, icell, bz2ibz, 

716 size, offset=(0, 0, 0), pad_width=2): 

717 """Interpolate values from Monkhorst-Pack sampling. 

718 

719 `monkhorst_pack_interpolate` takes a set of `values`, for example 

720 eigenvalues, that are resolved on a Monkhorst Pack k-point grid given by 

721 `size` and `offset` and interpolates the values onto the k-points 

722 in `path`. 

723 

724 Note 

725 ---- 

726 For the interpolation to work, path has to lie inside the domain 

727 that is spanned by the MP kpoint grid given by size and offset. 

728 

729 To try to ensure this we expand the domain slightly by adding additional 

730 entries along the edges and sides of the domain with values determined by 

731 wrapping the values to the opposite side of the domain. In this way we 

732 assume that the function to be interpolated is a periodic function in 

733 k-space. The padding width is determined by the `pad_width` parameter. 

734 

735 Parameters 

736 ---------- 

737 path: (nk, 3) array-like 

738 Desired path in units of reciprocal lattice vectors. 

739 values: (nibz, ...) array-like 

740 Values on Monkhorst-Pack grid. 

741 icell: (3, 3) array-like 

742 Reciprocal lattice vectors. 

743 bz2ibz: (nbz,) array-like of int 

744 Map from nbz points in BZ to nibz reduced points in IBZ. 

745 size: (3,) array-like of int 

746 Size of Monkhorst-Pack grid. 

747 offset: (3,) array-like 

748 Offset of Monkhorst-Pack grid. 

749 pad_width: int 

750 Padding width to aid interpolation 

751 

752 Returns 

753 ------- 

754 (nbz,) array-like 

755 *values* interpolated to *path*. 

756 """ 

757 from scipy.interpolate import LinearNDInterpolator 

758 

759 path = (np.asarray(path) + 0.5) % 1 - 0.5 

760 path = np.dot(path, icell) 

761 

762 # Fold out values from IBZ to BZ: 

763 v = np.asarray(values)[bz2ibz] 

764 v = v.reshape(tuple(size) + v.shape[1:]) 

765 

766 # Create padded Monkhorst-Pack grid: 

767 size = np.asarray(size) 

768 i = (np.indices(size + 2 * pad_width) 

769 .transpose((1, 2, 3, 0)).reshape((-1, 3))) 

770 k = (i - pad_width + 0.5) / size - 0.5 + offset 

771 k = np.dot(k, icell) 

772 

773 # Fill in boundary values: 

774 V = np.pad(v, [(pad_width, pad_width)] * 3 + 

775 [(0, 0)] * (v.ndim - 3), mode="wrap") 

776 

777 interpolate = LinearNDInterpolator(k, V.reshape((-1,) + V.shape[3:])) 

778 interpolated_points = interpolate(path) 

779 

780 # NaN values indicate points outside interpolation domain, if fail 

781 # try increasing padding 

782 assert not np.isnan(interpolated_points).any(), \ 

783 "Points outside interpolation domain. Try increasing pad_width." 

784 

785 return interpolated_points 

786 

787 

788# ChadiCohen k point grids. The k point grids are given in units of the 

789# reciprocal unit cell. The variables are named after the following 

790# convention: cc+'<Nkpoints>'+_+'shape'. For example an 18 k point 

791# sq(3)xsq(3) is named 'cc18_sq3xsq3'. 

792 

793cc6_1x1 = np.array([ 

794 1, 1, 0, 1, 0, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0, 

795 0, 1, 0]).reshape((6, 3)) / 3.0 

796 

797cc12_2x3 = np.array([ 

798 3, 4, 0, 3, 10, 0, 6, 8, 0, 3, -2, 0, 6, -4, 0, 

799 6, 2, 0, -3, 8, 0, -3, 2, 0, -3, -4, 0, -6, 4, 0, -6, -2, 0, -6, 

800 -8, 0]).reshape((12, 3)) / 18.0 

801 

802cc18_sq3xsq3 = np.array([ 

803 2, 2, 0, 4, 4, 0, 8, 2, 0, 4, -2, 0, 8, -4, 

804 0, 10, -2, 0, 10, -8, 0, 8, -10, 0, 2, -10, 0, 4, -8, 0, -2, -8, 

805 0, 2, -4, 0, -4, -4, 0, -2, -2, 0, -4, 2, 0, -2, 4, 0, -8, 4, 0, 

806 -4, 8, 0]).reshape((18, 3)) / 18.0 

807 

808cc18_1x1 = np.array([ 

809 2, 4, 0, 2, 10, 0, 4, 8, 0, 8, 4, 0, 8, 10, 0, 

810 10, 8, 0, 2, -2, 0, 4, -4, 0, 4, 2, 0, -2, 8, 0, -2, 2, 0, -2, -4, 

811 0, -4, 4, 0, -4, -2, 0, -4, -8, 0, -8, 2, 0, -8, -4, 0, -10, -2, 

812 0]).reshape((18, 3)) / 18.0 

813 

814cc54_sq3xsq3 = np.array([ 

815 4, -10, 0, 6, -10, 0, 0, -8, 0, 2, -8, 0, 6, 

816 -8, 0, 8, -8, 0, -4, -6, 0, -2, -6, 0, 2, -6, 0, 4, -6, 0, 8, -6, 

817 0, 10, -6, 0, -6, -4, 0, -2, -4, 0, 0, -4, 0, 4, -4, 0, 6, -4, 0, 

818 10, -4, 0, -6, -2, 0, -4, -2, 0, 0, -2, 0, 2, -2, 0, 6, -2, 0, 8, 

819 -2, 0, -8, 0, 0, -4, 0, 0, -2, 0, 0, 2, 0, 0, 4, 0, 0, 8, 0, 0, 

820 -8, 2, 0, -6, 2, 0, -2, 2, 0, 0, 2, 0, 4, 2, 0, 6, 2, 0, -10, 4, 

821 0, -6, 4, 0, -4, 4, 0, 0, 4, 0, 2, 4, 0, 6, 4, 0, -10, 6, 0, -8, 

822 6, 0, -4, 6, 0, -2, 6, 0, 2, 6, 0, 4, 6, 0, -8, 8, 0, -6, 8, 0, 

823 -2, 8, 0, 0, 8, 0, -6, 10, 0, -4, 10, 0]).reshape((54, 3)) / 18.0 

824 

825cc54_1x1 = np.array([ 

826 2, 2, 0, 4, 4, 0, 8, 8, 0, 6, 8, 0, 4, 6, 0, 6, 

827 10, 0, 4, 10, 0, 2, 6, 0, 2, 8, 0, 0, 2, 0, 0, 4, 0, 0, 8, 0, -2, 

828 6, 0, -2, 4, 0, -4, 6, 0, -6, 4, 0, -4, 2, 0, -6, 2, 0, -2, 0, 0, 

829 -4, 0, 0, -8, 0, 0, -8, -2, 0, -6, -2, 0, -10, -4, 0, -10, -6, 0, 

830 -6, -4, 0, -8, -6, 0, -2, -2, 0, -4, -4, 0, -8, -8, 0, 4, -2, 0, 

831 6, -2, 0, 6, -4, 0, 2, 0, 0, 4, 0, 0, 6, 2, 0, 6, 4, 0, 8, 6, 0, 

832 8, 0, 0, 8, 2, 0, 10, 4, 0, 10, 6, 0, 2, -4, 0, 2, -6, 0, 4, -6, 

833 0, 0, -2, 0, 0, -4, 0, -2, -6, 0, -4, -6, 0, -6, -8, 0, 0, -8, 0, 

834 -2, -8, 0, -4, -10, 0, -6, -10, 0]).reshape((54, 3)) / 18.0 

835 

836cc162_sq3xsq3 = np.array([ 

837 -8, 16, 0, -10, 14, 0, -7, 14, 0, -4, 14, 

838 0, -11, 13, 0, -8, 13, 0, -5, 13, 0, -2, 13, 0, -13, 11, 0, -10, 

839 11, 0, -7, 11, 0, -4, 11, 0, -1, 11, 0, 2, 11, 0, -14, 10, 0, -11, 

840 10, 0, -8, 10, 0, -5, 10, 0, -2, 10, 0, 1, 10, 0, 4, 10, 0, -16, 

841 8, 0, -13, 8, 0, -10, 8, 0, -7, 8, 0, -4, 8, 0, -1, 8, 0, 2, 8, 0, 

842 5, 8, 0, 8, 8, 0, -14, 7, 0, -11, 7, 0, -8, 7, 0, -5, 7, 0, -2, 7, 

843 0, 1, 7, 0, 4, 7, 0, 7, 7, 0, 10, 7, 0, -13, 5, 0, -10, 5, 0, -7, 

844 5, 0, -4, 5, 0, -1, 5, 0, 2, 5, 0, 5, 5, 0, 8, 5, 0, 11, 5, 0, 

845 -14, 4, 0, -11, 4, 0, -8, 4, 0, -5, 4, 0, -2, 4, 0, 1, 4, 0, 4, 4, 

846 0, 7, 4, 0, 10, 4, 0, -13, 2, 0, -10, 2, 0, -7, 2, 0, -4, 2, 0, 

847 -1, 2, 0, 2, 2, 0, 5, 2, 0, 8, 2, 0, 11, 2, 0, -11, 1, 0, -8, 1, 

848 0, -5, 1, 0, -2, 1, 0, 1, 1, 0, 4, 1, 0, 7, 1, 0, 10, 1, 0, 13, 1, 

849 0, -10, -1, 0, -7, -1, 0, -4, -1, 0, -1, -1, 0, 2, -1, 0, 5, -1, 

850 0, 8, -1, 0, 11, -1, 0, 14, -1, 0, -11, -2, 0, -8, -2, 0, -5, -2, 

851 0, -2, -2, 0, 1, -2, 0, 4, -2, 0, 7, -2, 0, 10, -2, 0, 13, -2, 0, 

852 -10, -4, 0, -7, -4, 0, -4, -4, 0, -1, -4, 0, 2, -4, 0, 5, -4, 0, 

853 8, -4, 0, 11, -4, 0, 14, -4, 0, -8, -5, 0, -5, -5, 0, -2, -5, 0, 

854 1, -5, 0, 4, -5, 0, 7, -5, 0, 10, -5, 0, 13, -5, 0, 16, -5, 0, -7, 

855 -7, 0, -4, -7, 0, -1, -7, 0, 2, -7, 0, 5, -7, 0, 8, -7, 0, 11, -7, 

856 0, 14, -7, 0, 17, -7, 0, -8, -8, 0, -5, -8, 0, -2, -8, 0, 1, -8, 

857 0, 4, -8, 0, 7, -8, 0, 10, -8, 0, 13, -8, 0, 16, -8, 0, -7, -10, 

858 0, -4, -10, 0, -1, -10, 0, 2, -10, 0, 5, -10, 0, 8, -10, 0, 11, 

859 -10, 0, 14, -10, 0, 17, -10, 0, -5, -11, 0, -2, -11, 0, 1, -11, 0, 

860 4, -11, 0, 7, -11, 0, 10, -11, 0, 13, -11, 0, 16, -11, 0, -1, -13, 

861 0, 2, -13, 0, 5, -13, 0, 8, -13, 0, 11, -13, 0, 14, -13, 0, 1, 

862 -14, 0, 4, -14, 0, 7, -14, 0, 10, -14, 0, 13, -14, 0, 5, -16, 0, 

863 8, -16, 0, 11, -16, 0, 7, -17, 0, 10, -17, 0]).reshape((162, 3)) / 27.0 

864 

865cc162_1x1 = np.array([ 

866 -8, -16, 0, -10, -14, 0, -7, -14, 0, -4, -14, 

867 0, -11, -13, 0, -8, -13, 0, -5, -13, 0, -2, -13, 0, -13, -11, 0, 

868 -10, -11, 0, -7, -11, 0, -4, -11, 0, -1, -11, 0, 2, -11, 0, -14, 

869 -10, 0, -11, -10, 0, -8, -10, 0, -5, -10, 0, -2, -10, 0, 1, -10, 

870 0, 4, -10, 0, -16, -8, 0, -13, -8, 0, -10, -8, 0, -7, -8, 0, -4, 

871 -8, 0, -1, -8, 0, 2, -8, 0, 5, -8, 0, 8, -8, 0, -14, -7, 0, -11, 

872 -7, 0, -8, -7, 0, -5, -7, 0, -2, -7, 0, 1, -7, 0, 4, -7, 0, 7, -7, 

873 0, 10, -7, 0, -13, -5, 0, -10, -5, 0, -7, -5, 0, -4, -5, 0, -1, 

874 -5, 0, 2, -5, 0, 5, -5, 0, 8, -5, 0, 11, -5, 0, -14, -4, 0, -11, 

875 -4, 0, -8, -4, 0, -5, -4, 0, -2, -4, 0, 1, -4, 0, 4, -4, 0, 7, -4, 

876 0, 10, -4, 0, -13, -2, 0, -10, -2, 0, -7, -2, 0, -4, -2, 0, -1, 

877 -2, 0, 2, -2, 0, 5, -2, 0, 8, -2, 0, 11, -2, 0, -11, -1, 0, -8, 

878 -1, 0, -5, -1, 0, -2, -1, 0, 1, -1, 0, 4, -1, 0, 7, -1, 0, 10, -1, 

879 0, 13, -1, 0, -10, 1, 0, -7, 1, 0, -4, 1, 0, -1, 1, 0, 2, 1, 0, 5, 

880 1, 0, 8, 1, 0, 11, 1, 0, 14, 1, 0, -11, 2, 0, -8, 2, 0, -5, 2, 0, 

881 -2, 2, 0, 1, 2, 0, 4, 2, 0, 7, 2, 0, 10, 2, 0, 13, 2, 0, -10, 4, 

882 0, -7, 4, 0, -4, 4, 0, -1, 4, 0, 2, 4, 0, 5, 4, 0, 8, 4, 0, 11, 4, 

883 0, 14, 4, 0, -8, 5, 0, -5, 5, 0, -2, 5, 0, 1, 5, 0, 4, 5, 0, 7, 5, 

884 0, 10, 5, 0, 13, 5, 0, 16, 5, 0, -7, 7, 0, -4, 7, 0, -1, 7, 0, 2, 

885 7, 0, 5, 7, 0, 8, 7, 0, 11, 7, 0, 14, 7, 0, 17, 7, 0, -8, 8, 0, 

886 -5, 8, 0, -2, 8, 0, 1, 8, 0, 4, 8, 0, 7, 8, 0, 10, 8, 0, 13, 8, 0, 

887 16, 8, 0, -7, 10, 0, -4, 10, 0, -1, 10, 0, 2, 10, 0, 5, 10, 0, 8, 

888 10, 0, 11, 10, 0, 14, 10, 0, 17, 10, 0, -5, 11, 0, -2, 11, 0, 1, 

889 11, 0, 4, 11, 0, 7, 11, 0, 10, 11, 0, 13, 11, 0, 16, 11, 0, -1, 

890 13, 0, 2, 13, 0, 5, 13, 0, 8, 13, 0, 11, 13, 0, 14, 13, 0, 1, 14, 

891 0, 4, 14, 0, 7, 14, 0, 10, 14, 0, 13, 14, 0, 5, 16, 0, 8, 16, 0, 

892 11, 16, 0, 7, 17, 0, 10, 17, 0]).reshape((162, 3)) / 27.0 

893 

894 

895# The following is a list of the critical points in the 1st Brillouin zone 

896# for some typical crystal structures following the conventions of Setyawan 

897# and Curtarolo [https://doi.org/10.1016/j.commatsci.2010.05.010]. 

898# 

899# In units of the reciprocal basis vectors. 

900# 

901# See http://en.wikipedia.org/wiki/Brillouin_zone 

902sc_special_points = { 

903 'cubic': {'G': [0, 0, 0], 

904 'M': [1 / 2, 1 / 2, 0], 

905 'R': [1 / 2, 1 / 2, 1 / 2], 

906 'X': [0, 1 / 2, 0]}, 

907 'fcc': {'G': [0, 0, 0], 

908 'K': [3 / 8, 3 / 8, 3 / 4], 

909 'L': [1 / 2, 1 / 2, 1 / 2], 

910 'U': [5 / 8, 1 / 4, 5 / 8], 

911 'W': [1 / 2, 1 / 4, 3 / 4], 

912 'X': [1 / 2, 0, 1 / 2]}, 

913 'bcc': {'G': [0, 0, 0], 

914 'H': [1 / 2, -1 / 2, 1 / 2], 

915 'P': [1 / 4, 1 / 4, 1 / 4], 

916 'N': [0, 0, 1 / 2]}, 

917 'tetragonal': {'G': [0, 0, 0], 

918 'A': [1 / 2, 1 / 2, 1 / 2], 

919 'M': [1 / 2, 1 / 2, 0], 

920 'R': [0, 1 / 2, 1 / 2], 

921 'X': [0, 1 / 2, 0], 

922 'Z': [0, 0, 1 / 2]}, 

923 'orthorhombic': {'G': [0, 0, 0], 

924 'R': [1 / 2, 1 / 2, 1 / 2], 

925 'S': [1 / 2, 1 / 2, 0], 

926 'T': [0, 1 / 2, 1 / 2], 

927 'U': [1 / 2, 0, 1 / 2], 

928 'X': [1 / 2, 0, 0], 

929 'Y': [0, 1 / 2, 0], 

930 'Z': [0, 0, 1 / 2]}, 

931 'hexagonal': {'G': [0, 0, 0], 

932 'A': [0, 0, 1 / 2], 

933 'H': [1 / 3, 1 / 3, 1 / 2], 

934 'K': [1 / 3, 1 / 3, 0], 

935 'L': [1 / 2, 0, 1 / 2], 

936 'M': [1 / 2, 0, 0]}} 

937 

938 

939# Old version of dictionary kept for backwards compatibility. 

940# Not for ordinary use. 

941ibz_points = {'cubic': {'Gamma': [0, 0, 0], 

942 'X': [0, 0 / 2, 1 / 2], 

943 'R': [1 / 2, 1 / 2, 1 / 2], 

944 'M': [0 / 2, 1 / 2, 1 / 2]}, 

945 'fcc': {'Gamma': [0, 0, 0], 

946 'X': [1 / 2, 0, 1 / 2], 

947 'W': [1 / 2, 1 / 4, 3 / 4], 

948 'K': [3 / 8, 3 / 8, 3 / 4], 

949 'U': [5 / 8, 1 / 4, 5 / 8], 

950 'L': [1 / 2, 1 / 2, 1 / 2]}, 

951 'bcc': {'Gamma': [0, 0, 0], 

952 'H': [1 / 2, -1 / 2, 1 / 2], 

953 'N': [0, 0, 1 / 2], 

954 'P': [1 / 4, 1 / 4, 1 / 4]}, 

955 'hexagonal': {'Gamma': [0, 0, 0], 

956 'M': [0, 1 / 2, 0], 

957 'K': [-1 / 3, 1 / 3, 0], 

958 'A': [0, 0, 1 / 2], 

959 'L': [0, 1 / 2, 1 / 2], 

960 'H': [-1 / 3, 1 / 3, 1 / 2]}, 

961 'tetragonal': {'Gamma': [0, 0, 0], 

962 'X': [1 / 2, 0, 0], 

963 'M': [1 / 2, 1 / 2, 0], 

964 'Z': [0, 0, 1 / 2], 

965 'R': [1 / 2, 0, 1 / 2], 

966 'A': [1 / 2, 1 / 2, 1 / 2]}, 

967 'orthorhombic': {'Gamma': [0, 0, 0], 

968 'R': [1 / 2, 1 / 2, 1 / 2], 

969 'S': [1 / 2, 1 / 2, 0], 

970 'T': [0, 1 / 2, 1 / 2], 

971 'U': [1 / 2, 0, 1 / 2], 

972 'X': [1 / 2, 0, 0], 

973 'Y': [0, 1 / 2, 0], 

974 'Z': [0, 0, 1 / 2]}}