Coverage for ase / dft / kpoints.py: 86.80%

1# fmt: off 

2from __future__ import annotations 

3 

4import re 

5import warnings 

6from typing import TYPE_CHECKING 

7 

8import numpy as np 

9 

10from ase import Atoms 

11from ase.cell import Cell 

12from ase.utils import jsonable 

13 

14if TYPE_CHECKING: 

15 from ase.spectrum.band_structure import BandStructure 

16 

17 

18def monkhorst_pack(size): 

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

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

21 raise ValueError(f'Illegal size: {list(size)}') 

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

23 return (kpts + 0.5) / size - 0.5 

24 

25 

26def get_monkhorst_pack_size_and_offset(kpts): 

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

28 

29 Returns (size, offset), where:: 

30 

31 kpts = monkhorst_pack(size) + offset. 

32 

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

34 

35 if len(kpts) == 1: 

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

37 

38 size = np.zeros(3, int) 

39 for c in range(3): 

40 # Determine increment between k-points along current axis 

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

42 

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

44 if delta > 1e-8: 

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

46 else: 

47 size[c] = 1 

48 

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

50 kpts0 = monkhorst_pack(size) 

51 offsets = kpts - kpts0 

52 

53 # All offsets must be identical: 

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

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

56 

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

58 

59 

60def mindistance2monkhorstpack(atoms, *, 

61 min_distance, 

62 maxperdim=16, 

63 even=True): 

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

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

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

67 (nx, ny, nz)-supercell. 

68 

69 Compared to ase.calculators.calculator kptdensity2monkhorstpack 

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

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

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

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

74 invariant to particular choice of cell representations. 

75 

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

77 """ 

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

79 min_distance, maxperdim, even) 

80 

81 

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

83 from ase.neighborlist import NeighborList 

84 

85 step = 2 if even else 1 

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

87 self_interaction=False, bothways=False) 

88 

89 def check(nkpts_c): 

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

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

92 

93 def generate_mpgrids(): 

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

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

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

97 yield from sorted(nkpts_nc, key=np.prod) 

98 

99 try: 

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

101 except StopIteration: 

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

103 ' Try running with a larger maxperdim.') 

104 

105 

106def get_monkhorst_shape(kpts): 

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

108 return get_monkhorst_pack_size_and_offset(kpts)[0] 

109 

110 

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

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

113 

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

115 coordinates, the other is determined. 

116 

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

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

119 """ 

120 if ckpts_kv is None: 

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

122 return np.dot(skpts_kc, icell_cv) 

123 elif skpts_kc is None: 

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

125 else: 

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

127 

128 

129def parse_path_string(s): 

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

131 

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

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

134 list of labels. 

135 

136 Examples 

137 -------- 

138 

139 >>> parse_path_string('GX') 

140 [['G', 'X']] 

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

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

143 """ 

144 paths = [] 

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

146 names = [name 

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

148 if name] 

149 paths.append(names) 

150 return paths 

151 

152 

153def resolve_kpt_path_string(path, special_points): 

154 paths = parse_path_string(path) 

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

156 for subpath in paths] 

157 return paths, coords 

158 

159 

160def resolve_custom_points(pathspec, special_points, eps): 

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

162 

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

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

165 

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

167 are generated dynamically as necessary, updating the special_point 

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

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

170 labelled as such.""" 

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

172 # be lazy and call it on scaled ones. 

173 

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

175 special_points = dict(special_points) 

176 

177 if len(pathspec) == 0: 

178 return '', special_points 

179 

180 if isinstance(pathspec, str): 

181 pathspec = parse_path_string(pathspec) 

182 

183 def looks_like_single_kpoint(obj): 

184 if isinstance(obj, str): 

185 return True 

186 try: 

187 arr = np.asarray(obj, float) 

188 except ValueError: 

189 return False 

190 else: 

191 return arr.shape == (3,) 

192 

193 # We accept inputs that are either 

194 # 1) string notation 

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

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

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

198 for thing in pathspec: 

199 if looks_like_single_kpoint(thing): 

200 pathspec = [pathspec] 

201 break 

202 

203 def name_generator(): 

204 counter = 0 

205 while True: 

206 name = f'Kpt{counter}' 

207 yield name 

208 counter += 1 

209 

210 custom_names = name_generator() 

211 

212 labelseq = [] 

213 for subpath in pathspec: 

214 for kpt in subpath: 

215 if isinstance(kpt, str): 

216 if kpt not in special_points: 

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

218 .format(kpt, 

219 ''.join(special_points))) 

220 labelseq.append(kpt) 

221 continue 

222 

223 kpt = np.asarray(kpt, float) 

224 if kpt.shape != (3,): 

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

226 

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

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

229 labelseq.append(key) 

230 break # Already present 

231 else: 

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

233 name = next(custom_names) 

234 while name in special_points: 

235 name = next(custom_names) 

236 special_points[name] = kpt 

237 labelseq.append(name) 

238 labelseq.append(',') 

239 

240 last = labelseq.pop() 

241 assert last == ',' 

242 return ''.join(labelseq), special_points 

243 

244 

245def normalize_special_points(special_points): 

246 dct = {} 

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

248 if not isinstance(name, str): 

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

250 if np.shape(value) != (3,): 

251 raise ValueError('Expected 3 kpoint coordinates') 

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

253 return dct 

254 

255 

256@jsonable('bandpath') 

257class BandPath: 

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

259 

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

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

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

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

264 

265 >>> from ase.lattice import CUB 

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

267 >>> path 

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

269 

270 Band paths support JSON I/O: 

271 

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

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

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

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

276 

277 """ 

278 

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

280 special_points=None, path=None): 

281 if kpts is None: 

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

283 

284 if special_points is None: 

285 special_points = {} 

286 else: 

287 special_points = normalize_special_points(special_points) 

288 

289 if path is None: 

290 path = '' 

291 

292 cell = Cell(cell) 

293 self._cell = cell 

294 kpts = np.asarray(kpts) 

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

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

297 self._kpts = kpts 

298 self._special_points = special_points 

299 if not isinstance(path, str): 

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

301 self._path = path 

302 

303 @property 

304 def cell(self) -> Cell: 

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

306 return self._cell 

307 

308 @property 

309 def icell(self) -> Cell: 

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

311 return self._icell 

312 

313 @property 

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

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

316 

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

318 return self._kpts 

319 

320 @property 

321 def special_points(self) -> dict[str, np.ndarray]: 

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

323 

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

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

326 

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

328 consider deepcopying it.""" 

329 return self._special_points 

330 

331 @property 

332 def path(self) -> str: 

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

334 

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

336 

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

338 return self._path 

339 

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

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

342 

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

344 The operation will typically be a permutation/flipping 

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

346 # XXX acceptable operations are probably only those 

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

348 # 

349 # We should insert a check. 

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

351 # if they change lengths, volume etc. 

352 special_points = { 

353 name: value @ op for name, value in self.special_points.items() 

354 } 

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

356 special_points=special_points, 

357 path=self.path) 

358 

359 def todict(self): 

360 return {'kpts': self.kpts, 

361 'special_points': self.special_points, 

362 'labelseq': self.path, 

363 'cell': self.cell} 

364 

365 def interpolate( 

366 self, 

367 path: str | None = None, 

368 npoints: int | None = None, 

369 special_points: dict[str, np.ndarray] | None = None, 

370 density: float | None = None, 

371 ) -> BandPath: 

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

373 if path is None: 

374 path = self.path 

375 

376 if special_points is None: 

377 special_points = self.special_points 

378 

379 _pathnames, pathcoords = resolve_kpt_path_string(path, special_points) 

380 kpts, _x, _X = paths2kpts(pathcoords, self.cell, npoints, density) 

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

382 special_points=special_points) 

383 

384 def _scale(self, coords): 

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

386 

387 def __repr__(self): 

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

389 .format(self.__class__.__name__, 

390 repr(self.path), 

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

392 len(self.kpts))) 

393 

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

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

396 return self._scale(self.kpts) 

397 

398 def __iter__(self): 

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

400 

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

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

403 special_x_coords), and people would use tuple unpacking for 

404 those. 

405 

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

407 It will be removed in the future. 

408 

409 """ 

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

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

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

413 yield self.kpts 

414 

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

416 yield x 

417 yield xspecial 

418 

419 def __getitem__(self, index): 

420 # Temp compatibility stuff, see __iter__ 

421 return tuple(self)[index] 

422 

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

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

425 

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

427 

428 index2name = self._find_special_point_indices(eps) 

429 indices = sorted(index2name) 

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

431 xcoords, special_xcoords = indices_to_axis_coords( 

432 indices, self.kpts, self.cell) 

433 return xcoords, special_xcoords, labels 

434 

435 def _find_special_point_indices(self, eps): 

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

437 

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

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

440 # to fully make sense! 

441 index2name = {} 

442 nkpts = len(self.kpts) 

443 

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

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

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

447 args = np.argwhere(distances < eps) 

448 for arg in args.flat: 

449 dist = distances[arg] 

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

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

452 if not any(neighbours < dist): 

453 index2name[arg] = name 

454 

455 return index2name 

456 

457 def plot(self, **plotkwargs): 

458 """Visualize this bandpath. 

459 

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

461 import ase.dft.bz as bz 

462 

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

464 plotkwargs.pop('dimension', None) 

465 

466 special_points = self.special_points 

467 labelseq, coords = resolve_kpt_path_string(self.path, 

468 special_points) 

469 

470 paths = [] 

471 points_already_plotted = set() 

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

473 subpath_coords = np.array(subpath_coords) 

474 points_already_plotted.update(subpath_labels) 

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

476 

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

478 # not plotted already: 

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

480 if label not in points_already_plotted: 

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

482 

483 kw = {'vectors': True, 

484 'pointstyle': {'marker': '.'}} 

485 

486 kw.update(plotkwargs) 

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

488 points=self.cartesian_kpts(), 

489 **kw) 

490 

491 def free_electron_band_structure(self, **kwargs) -> BandStructure: 

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

493 

494 Keyword arguments are passed to 

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

496 

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

498 from ase.calculators.test import FreeElectrons 

499 from ase.spectrum.band_structure import calculate_band_structure 

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

501 atoms.calc = FreeElectrons(**kwargs) 

502 bs = calculate_band_structure(atoms, path=self) 

503 return bs 

504 

505 

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

507 eps=2e-4): 

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

509 

510 path: list or str 

511 Can be: 

512 

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

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

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

516 cell: 3x3 

517 Unit cell of the atoms. 

518 npoints: int 

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

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

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

522 density: float 

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

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

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

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

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

528 will be used. 

529 special_points: dict or None 

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

531 points will be derived from the cell. 

532 eps: float 

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

534 

535 You may define npoints or density but not both. 

536 

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

538 

539 cell = Cell.ascell(cell) 

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

541 special_points=special_points, eps=eps) 

542 

543 

544DEFAULT_KPTS_DENSITY = 5 # points per 1/Angstrom 

545 

546 

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

548 if npoints is not None and density is not None: 

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

550 points = np.concatenate(paths) 

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

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

553 

554 i = 0 

555 for path in paths[:-1]: 

556 i += len(path) 

557 lengths[i - 1] = 0 

558 

559 length = sum(lengths) 

560 

561 if npoints is None: 

562 if density is None: 

563 density = DEFAULT_KPTS_DENSITY 

564 # Set npoints using the length of the path 

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

566 

567 kpts = [] 

568 x0 = 0 

569 x = [] 

570 X = [0] 

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

572 diff = length - x0 

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

574 n = 0 

575 else: 

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

577 

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

579 kpts.append(P + t * d) 

580 x.append(x0 + t * L) 

581 x0 += L 

582 X.append(x0) 

583 if len(points): 

584 kpts.append(points[-1]) 

585 x.append(x0) 

586 

587 if len(kpts) == 0: 

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

589 

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

591 

592 

593get_bandpath = bandpath # old name 

594 

595 

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

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

598 # XXX Should use the Cartesian kpoints. 

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

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

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

602 N = len(kpts) 

603 indices = [] 

604 if N > 0: 

605 indices.append(0) 

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

607 indices.append(N - 1) 

608 return indices 

609 

610 

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

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

613 

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

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

616 as xticklabels. 

617 

618 Parameters 

619 ---------- 

620 

621 kpts: list 

622 List of scaled k-points. 

623 

624 cell: list 

625 Unit cell of the atomic structure. 

626 

627 Returns 

628 ------- 

629 

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

631 the second is x coordinates of the special points, 

632 the third is the special points as strings. 

633 """ 

634 

635 if special_points is None: 

636 special_points = get_special_points(cell) 

637 points = np.asarray(kpts) 

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

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

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

641 

642 labels = [] 

643 for kpt in points[indices]: 

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

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

646 break 

647 else: 

648 # No exact match. Try modulus 1: 

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

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

651 break 

652 else: 

653 label = '?' 

654 labels.append(label) 

655 

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

657 return xcoords, ixcoords, labels 

658 

659 

660def indices_to_axis_coords(indices, points, cell): 

661 jump = False # marks a discontinuity in the path 

662 xcoords = [0] 

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

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

665 length = 0 

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

667 else: 

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

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

670 jump = False 

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

672 

673 xcoords = np.array(xcoords) 

674 return xcoords, xcoords[indices] 

675 

676 

677special_paths = { 

678 'cubic': 'GXMGRX,MR', 

679 'fcc': 'GXWKGLUWLK,UX', 

680 'bcc': 'GHNGPH,PN', 

681 'tetragonal': 'GXMGZRAZXR,MA', 

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

683 'hexagonal': 'GMKGALHA,LM,KH', 

684 'monoclinic': 'GYHCEM1AXH1,MDZ,YD', 

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

686 'rhombohedral type 2': 'GPZQGFP1Q1LZ'} 

687 

688 

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

690 """Return dict of special points. 

691 

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

693 Curtarolo:: 

694 

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

696 

697 cell: 3x3 ndarray 

698 Unit cell. 

699 lattice: str 

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

701 bcc, orthorhombic, tetragonal, hexagonal or monoclinic. 

702 eps: float 

703 Tolerance for cell-check. 

704 """ 

705 

706 if isinstance(cell, str): 

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

708 'argument') 

709 lattice, cell = cell, lattice 

710 

711 cell = Cell.ascell(cell) 

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

713 # from the canonical cell to the given one. 

714 # 

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

716 path = cell.bandpath(npoints=0) 

717 return path.special_points 

718 

719 

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

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

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

723 

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

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

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

727 in `path`. 

728 

729 Note 

730 ---- 

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

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

733 

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

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

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

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

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

739 

740 Parameters 

741 ---------- 

742 path: (nk, 3) array-like 

743 Desired path in units of reciprocal lattice vectors. 

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

745 Values on Monkhorst-Pack grid. 

746 icell: (3, 3) array-like 

747 Reciprocal lattice vectors. 

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

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

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

751 Size of Monkhorst-Pack grid. 

752 offset: (3,) array-like 

753 Offset of Monkhorst-Pack grid. 

754 pad_width: int 

755 Padding width to aid interpolation 

756 

757 Returns 

758 ------- 

759 (nbz,) array-like 

760 *values* interpolated to *path*. 

761 """ 

762 from scipy.interpolate import LinearNDInterpolator 

763 

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

765 path = np.dot(path, icell) 

766 

767 # Fold out values from IBZ to BZ: 

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

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

770 

771 # Create padded Monkhorst-Pack grid: 

772 size = np.asarray(size) 

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

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

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

776 k = np.dot(k, icell) 

777 

778 # Fill in boundary values: 

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

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

781 

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

783 interpolated_points = interpolate(path) 

784 

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

786 # try increasing padding 

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

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

789 

790 return interpolated_points 

791 

792 

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

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

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

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

797 

798cc6_1x1 = np.array([ 

799 1, 1, 0, 1, 0, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0, 

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

801 

802cc12_2x3 = np.array([ 

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

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

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

806 

807cc18_sq3xsq3 = np.array([ 

808 2, 2, 0, 4, 4, 0, 8, 2, 0, 4, -2, 0, 8, -4, 

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

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

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

812 

813cc18_1x1 = np.array([ 

814 2, 4, 0, 2, 10, 0, 4, 8, 0, 8, 4, 0, 8, 10, 0, 

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

816 0, -4, 4, 0, -4, -2, 0, -4, -8, 0, -8, 2, 0, -8, -4, 0, -10, -2, 

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

818 

819cc54_sq3xsq3 = np.array([ 

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

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