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 >>> parse_path_string('GX') 

139 [['G', 'X']] 

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

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

142 """ 

143 paths = [] 

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

145 names = [name 

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

147 if name] 

148 paths.append(names) 

149 return paths 

150 

151 

152def resolve_kpt_path_string(path, special_points): 

153 paths = parse_path_string(path) 

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

155 for subpath in paths] 

156 return paths, coords 

157 

158 

159def resolve_custom_points(pathspec, special_points, eps): 

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

161 

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

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

164 

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

166 are generated dynamically as necessary, updating the special_point 

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

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

169 labelled as such.""" 

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

171 # be lazy and call it on scaled ones. 

172 

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

174 special_points = dict(special_points) 

175 

176 if len(pathspec) == 0: 

177 return '', special_points 

178 

179 if isinstance(pathspec, str): 

180 pathspec = parse_path_string(pathspec) 

181 

182 def looks_like_single_kpoint(obj): 

183 if isinstance(obj, str): 

184 return True 

185 try: 

186 arr = np.asarray(obj, float) 

187 except ValueError: 

188 return False 

189 else: 

190 return arr.shape == (3,) 

191 

192 # We accept inputs that are either 

193 # 1) string notation 

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

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

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

197 for thing in pathspec: 

198 if looks_like_single_kpoint(thing): 

199 pathspec = [pathspec] 

200 break 

201 

202 def name_generator(): 

203 counter = 0 

204 while True: 

205 name = f'Kpt{counter}' 

206 yield name 

207 counter += 1 

208 

209 custom_names = name_generator() 

210 

211 labelseq = [] 

212 for subpath in pathspec: 

213 for kpt in subpath: 

214 if isinstance(kpt, str): 

215 if kpt not in special_points: 

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

217 .format(kpt, 

218 ''.join(special_points))) 

219 labelseq.append(kpt) 

220 continue 

221 

222 kpt = np.asarray(kpt, float) 

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

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

225 

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

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

228 labelseq.append(key) 

229 break # Already present 

230 else: 

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

232 name = next(custom_names) 

233 while name in special_points: 

234 name = next(custom_names) 

235 special_points[name] = kpt 

236 labelseq.append(name) 

237 labelseq.append(',') 

238 

239 last = labelseq.pop() 

240 assert last == ',' 

241 return ''.join(labelseq), special_points 

242 

243 

244def normalize_special_points(special_points): 

245 dct = {} 

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

247 if not isinstance(name, str): 

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

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

250 raise ValueError('Expected 3 kpoint coordinates') 

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

252 return dct 

253 

254 

255@jsonable('bandpath') 

256class BandPath: 

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

258 

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

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

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

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

263 

264 >>> from ase.lattice import CUB 

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

266 >>> path 

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

268 

269 Band paths support JSON I/O: 

270 

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

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

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

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

275 

276 """ 

277 

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

279 special_points=None, path=None): 

280 if kpts is None: 

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

282 

283 if special_points is None: 

284 special_points = {} 

285 else: 

286 special_points = normalize_special_points(special_points) 

287 

288 if path is None: 

289 path = '' 

290 

291 cell = Cell(cell) 

292 self._cell = cell 

293 kpts = np.asarray(kpts) 

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

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

296 self._kpts = kpts 

297 self._special_points = special_points 

298 if not isinstance(path, str): 

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

300 self._path = path 

301 

302 @property 

303 def cell(self) -> Cell: 

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

305 return self._cell 

306 

307 @property 

308 def icell(self) -> Cell: 

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

310 return self._icell 

311 

312 @property 

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

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

315 

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

317 return self._kpts 

318 

319 @property 

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

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

322 

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

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

325 

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

327 consider deepcopying it.""" 

328 return self._special_points 

329 

330 @property 

331 def path(self) -> str: 

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

333 

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

335 

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

337 return self._path 

338 

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

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

341 

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

343 The operation will typically be a permutation/flipping 

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

345 # XXX acceptable operations are probably only those 

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

347 # 

348 # We should insert a check. 

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

350 # if they change lengths, volume etc. 

351 special_points = { 

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

353 } 

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

355 special_points=special_points, 

356 path=self.path) 

357 

358 def todict(self): 

359 return {'kpts': self.kpts, 

360 'special_points': self.special_points, 

361 'labelseq': self.path, 

362 'cell': self.cell} 

363 

364 def interpolate( 

365 self, 

366 path: str = None, 

367 npoints: int = None, 

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

369 density: float = None, 

370 ) -> BandPath: 

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

372 if path is None: 

373 path = self.path 

374 

375 if special_points is None: 

376 special_points = self.special_points 

377 

378 _pathnames, pathcoords = resolve_kpt_path_string(path, special_points) 

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

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

381 special_points=special_points) 

382 

383 def _scale(self, coords): 

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

385 

386 def __repr__(self): 

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

388 .format(self.__class__.__name__, 

389 repr(self.path), 

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

391 len(self.kpts))) 

392 

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

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

395 return self._scale(self.kpts) 

396 

397 def __iter__(self): 

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

399 

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

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

402 special_x_coords), and people would use tuple unpacking for 

403 those. 

404 

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

406 It will be removed in the future. 

407 

408 """ 

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

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

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

412 yield self.kpts 

413 

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

415 yield x 

416 yield xspecial 

417 

418 def __getitem__(self, index): 

419 # Temp compatibility stuff, see __iter__ 

420 return tuple(self)[index] 

421 

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

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

424 

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

426 

427 index2name = self._find_special_point_indices(eps) 

428 indices = sorted(index2name) 

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

430 xcoords, special_xcoords = indices_to_axis_coords( 

431 indices, self.kpts, self.cell) 

432 return xcoords, special_xcoords, labels 

433 

434 def _find_special_point_indices(self, eps): 

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

436 

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

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

439 # to fully make sense! 

440 index2name = {} 

441 nkpts = len(self.kpts) 

442 

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

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

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

446 args = np.argwhere(distances < eps) 

447 for arg in args.flat: 

448 dist = distances[arg] 

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

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

451 if not any(neighbours < dist): 

452 index2name[arg] = name 

453 

454 return index2name 

455 

456 def plot(self, **plotkwargs): 

457 """Visualize this bandpath. 

458 

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

460 import ase.dft.bz as bz 

461 

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

463 plotkwargs.pop('dimension', None) 

464 

465 special_points = self.special_points 

466 labelseq, coords = resolve_kpt_path_string(self.path, 

467 special_points) 

468 

469 paths = [] 

470 points_already_plotted = set() 

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

472 subpath_coords = np.array(subpath_coords) 

473 points_already_plotted.update(subpath_labels) 

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

475 

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

477 # not plotted already: 

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

479 if label not in points_already_plotted: 

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

481 

482 kw = {'vectors': True, 

483 'pointstyle': {'marker': '.'}} 

484 

485 kw.update(plotkwargs) 

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

487 points=self.cartesian_kpts(), 

488 **kw) 

489 

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

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

492 

493 Keyword arguments are passed to 

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

495 

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

497 from ase.calculators.test import FreeElectrons 

498 from ase.spectrum.band_structure import calculate_band_structure 

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

500 atoms.calc = FreeElectrons(**kwargs) 

501 bs = calculate_band_structure(atoms, path=self) 

502 return bs 

503 

504 

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

506 eps=2e-4): 

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

508 

509 path: list or str 

510 Can be: 

511 

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

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

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

515 cell: 3x3 

516 Unit cell of the atoms. 

517 npoints: int 

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

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

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

521 density: float 

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

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

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

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

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

527 will be used. 

528 special_points: dict or None 

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

530 points will be derived from the cell. 

531 eps: float 

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

533 

534 You may define npoints or density but not both. 

535 

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

537 

538 cell = Cell.ascell(cell) 

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

540 special_points=special_points, eps=eps) 

541 

542 

543DEFAULT_KPTS_DENSITY = 5 # points per 1/Angstrom 

544 

545 

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

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

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

549 points = np.concatenate(paths) 

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

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

552 

553 i = 0 

554 for path in paths[:-1]: 

555 i += len(path) 

556 lengths[i - 1] = 0 

557 

558 length = sum(lengths) 

559 

560 if npoints is None: 

561 if density is None: 

562 density = DEFAULT_KPTS_DENSITY 

563 # Set npoints using the length of the path 

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

565 

566 kpts = [] 

567 x0 = 0 

568 x = [] 

569 X = [0] 

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

571 diff = length - x0 

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

573 n = 0 

574 else: 

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

576 

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

578 kpts.append(P + t * d) 

579 x.append(x0 + t * L) 

580 x0 += L 

581 X.append(x0) 

582 if len(points): 

583 kpts.append(points[-1]) 

584 x.append(x0) 

585 

586 if len(kpts) == 0: 

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

588 

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

590 

591 

592get_bandpath = bandpath # old name 

593 

594 

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

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

597 # XXX Should use the Cartesian kpoints. 

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

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

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

601 N = len(kpts) 

602 indices = [] 

603 if N > 0: 

604 indices.append(0) 

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

606 indices.append(N - 1) 

607 return indices 

608 

609 

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

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

612 

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

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

615 as xticklabels. 

616 

617 Parameters: 

618 

619 kpts: list 

620 List of scaled k-points. 

621 

622 cell: list 

623 Unit cell of the atomic structure. 

624 

625 Returns: 

626 

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

628 the second is x coordinates of the special points, 

629 the third is the special points as strings. 

630 """ 

631 

632 if special_points is None: 

633 special_points = get_special_points(cell) 

634 points = np.asarray(kpts) 

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

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

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

638 

639 labels = [] 

640 for kpt in points[indices]: 

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

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

643 break 

644 else: 

645 # No exact match. Try modulus 1: 

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

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

648 break 

649 else: 

650 label = '?' 

651 labels.append(label) 

652 

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

654 return xcoords, ixcoords, labels 

655 

656 

657def indices_to_axis_coords(indices, points, cell): 

658 jump = False # marks a discontinuity in the path 

659 xcoords = [0] 

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

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

662 length = 0 

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

664 else: 

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

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

667 jump = False 

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

669 

670 xcoords = np.array(xcoords) 

671 return xcoords, xcoords[indices] 

672 

673 

674special_paths = { 

675 'cubic': 'GXMGRX,MR', 

676 'fcc': 'GXWKGLUWLK,UX', 

677 'bcc': 'GHNGPH,PN', 

678 'tetragonal': 'GXMGZRAZXR,MA', 

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

680 'hexagonal': 'GMKGALHA,LM,KH', 

681 'monoclinic': 'GYHCEM1AXH1,MDZ,YD', 

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

683 'rhombohedral type 2': 'GPZQGFP1Q1LZ'} 

684 

685 

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

687 """Return dict of special points. 

688 

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

690 Curtarolo:: 

691 

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

693 

694 cell: 3x3 ndarray 

695 Unit cell. 

696 lattice: str 

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

698 bcc, orthorhombic, tetragonal, hexagonal or monoclinic. 

699 eps: float 

700 Tolerance for cell-check. 

701 """ 

702 

703 if isinstance(cell, str): 

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

705 'argument') 

706 lattice, cell = cell, lattice 

707 

708 cell = Cell.ascell(cell) 

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

710 # from the canonical cell to the given one. 

711 # 

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

713 path = cell.bandpath(npoints=0) 

714 return path.special_points 

715 

716 

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

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

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

720 

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

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

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

724 in `path`. 

725 

726 Note 

727 ---- 

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

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

730 

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

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

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

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

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

736 

737 Parameters 

738 ---------- 

739 path: (nk, 3) array-like 

740 Desired path in units of reciprocal lattice vectors. 

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

742 Values on Monkhorst-Pack grid. 

743 icell: (3, 3) array-like 

744 Reciprocal lattice vectors. 

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

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

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

748 Size of Monkhorst-Pack grid. 

749 offset: (3,) array-like 

750 Offset of Monkhorst-Pack grid. 

751 pad_width: int 

752 Padding width to aid interpolation 

753 

754 Returns 

755 ------- 

756 (nbz,) array-like 

757 *values* interpolated to *path*. 

758 """ 

759 from scipy.interpolate import LinearNDInterpolator 

760 

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

762 path = np.dot(path, icell) 

763 

764 # Fold out values from IBZ to BZ: 

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

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

767 

768 # Create padded Monkhorst-Pack grid: 

769 size = np.asarray(size) 

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

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

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

773 k = np.dot(k, icell) 

774 

775 # Fill in boundary values: 

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

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

778 

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

780 interpolated_points = interpolate(path) 

781 

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

783 # try increasing padding 

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

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

786 

787 return interpolated_points 

788 

789 

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

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

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

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

794 

795cc6_1x1 = np.array([ 

796 1, 1, 0, 1, 0, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0, 

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

798 

799cc12_2x3 = np.array([ 

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

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

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

803 

804cc18_sq3xsq3 = np.array([ 

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

806 0, 10, -2, 0, 10, -8, 0, 8, -10, 0, 2, -10, 0, 4, -8, 0, -2, -8, 

807 0, 2, -4, 0, -4, -4, 0, -2, -2, 0, -4, 2, 0, -2, 4, 0, -8, 4, 0, 

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

809 

810cc18_1x1 = np.array([ 

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

812 10, 8, 0, 2, -2, 0, 4, -4, 0, 4, 2, 0, -2, 8, 0, -2, 2, 0, -2, -4, 

813 0, -4, 4, 0, -4, -2, 0, -4, -8, 0, -8, 2, 0, -8, -4, 0, -10, -2, 

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

815 

816cc54_sq3xsq3 = np.array([ 

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

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

819 0, 10, -6, 0, -6, -4, 0, -2, -4, 0, 0, -4, 0, 4, -4, 0, 6, -4, 0,