Coverage for /builds/kinetik161/ase/ase/dft/kpoints.py: 86.67%

1import re 

2import warnings 

3from typing import Dict 

4 

5import numpy as np 

6 

7import ase # Annotations 

8from ase.cell import Cell 

9from ase.utils import jsonable 

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(f'Illegal size: {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=np.prod) 

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 = f'Kpt{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 

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

273 special_points=None, path=None): 

274 if kpts is None: 

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

276 

277 if special_points is None: 

278 special_points = {} 

279 else: 

280 special_points = normalize_special_points(special_points) 

281 

282 if path is None: 

283 path = '' 

284 

285 cell = Cell(cell) 

286 self._cell = cell 

287 kpts = np.asarray(kpts) 

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

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

290 self._kpts = kpts 

291 self._special_points = special_points 

292 if not isinstance(path, str): 

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

294 self._path = path 

295 

296 @property 

297 def cell(self) -> Cell: 

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

299 return self._cell 

300 

301 @property 

302 def icell(self) -> Cell: 

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

304 return self._icell 

305 

306 @property 

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

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

309 

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

311 return self._kpts 

312 

313 @property 

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

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

316 

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

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

319 

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

321 consider deepcopying it.""" 

322 return self._special_points 

323 

324 @property 

325 def path(self) -> str: 

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

327 

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

329 

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

331 return self._path 

332 

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

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

335 

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

337 The operation will typically be a permutation/flipping 

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

339 # XXX acceptable operations are probably only those 

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

341 # 

342 # We should insert a check. 

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

344 # if they change lengths, volume etc. 

345 special_points = {} 

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

347 special_points[name] = value @ op 

348 

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

350 special_points=special_points, 

351 path=self.path) 

352 

353 def todict(self): 

354 return {'kpts': self.kpts, 

355 'special_points': self.special_points, 

356 'labelseq': self.path, 

357 'cell': self.cell} 

358 

359 def interpolate( 

360 self, 

361 path: str = None, 

362 npoints: int = None, 

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

364 density: float = None, 

365 ) -> 'BandPath': 

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

367 if path is None: 

368 path = self.path 

369 

370 if special_points is None: 

371 special_points = self.special_points 

372 

373 pathnames, pathcoords = resolve_kpt_path_string(path, special_points) 

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

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

376 special_points=special_points) 

377 

378 def _scale(self, coords): 

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

380 

381 def __repr__(self): 

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

383 .format(self.__class__.__name__, 

384 repr(self.path), 

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

386 len(self.kpts))) 

387 

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

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

390 return self._scale(self.kpts) 

391 

392 def __iter__(self): 

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

394 

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

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

397 special_x_coords), and people would use tuple unpacking for 

398 those. 

399 

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

401 It will be removed in the future. 

402 

403 """ 

404 import warnings 

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

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

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

408 yield self.kpts 

409 

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

411 yield x 

412 yield xspecial 

413 

414 def __getitem__(self, index): 

415 # Temp compatibility stuff, see __iter__ 

416 return tuple(self)[index] 

417 

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

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

420 

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

422 

423 index2name = self._find_special_point_indices(eps) 

424 indices = sorted(index2name) 

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

426 xcoords, special_xcoords = indices_to_axis_coords( 

427 indices, self.kpts, self.cell) 

428 return xcoords, special_xcoords, labels 

429 

430 def _find_special_point_indices(self, eps): 

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

432 

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

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

435 # to fully make sense! 

436 index2name = {} 

437 nkpts = len(self.kpts) 

438 

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

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

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

442 args = np.argwhere(distances < eps) 

443 for arg in args.flat: 

444 dist = distances[arg] 

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

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

447 if not any(neighbours < dist): 

448 index2name[arg] = name 

449 

450 return index2name 

451 

452 def plot(self, **plotkwargs): 

453 """Visualize this bandpath. 

454 

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

456 import ase.dft.bz as bz 

457 

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

459 plotkwargs.pop('dimension', None) 

460 

461 special_points = self.special_points 

462 labelseq, coords = resolve_kpt_path_string(self.path, 

463 special_points) 

464 

465 paths = [] 

466 points_already_plotted = set() 

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

468 subpath_coords = np.array(subpath_coords) 

469 points_already_plotted.update(subpath_labels) 

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

471 

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

473 # not plotted already: 

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

475 if label not in points_already_plotted: 

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

477 

478 kw = {'vectors': True, 

479 'pointstyle': {'marker': '.'}} 

480 

481 kw.update(plotkwargs) 

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

483 points=self.cartesian_kpts(), 

484 **kw) 

485 

486 def free_electron_band_structure( 

487 self, **kwargs 

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

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

490 

491 Keyword arguments are passed to 

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

493 

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

495 from ase import Atoms 

496 from ase.calculators.test import FreeElectrons 

497 from ase.spectrum.band_structure import calculate_band_structure 

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

499 atoms.calc = FreeElectrons(**kwargs) 

500 bs = calculate_band_structure(atoms, path=self) 

501 return bs 

502 

503 

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

505 eps=2e-4): 

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

507 

508 path: list or str 

509 Can be: 

510 

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

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

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

514 cell: 3x3 

515 Unit cell of the atoms. 

516 npoints: int 

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

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

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

520 density: float 

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

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

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

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

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

526 will be used. 

527 special_points: dict or None 

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

529 points will be derived from the cell. 

530 eps: float 

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

532 

533 You may define npoints or density but not both. 

534 

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

536 

537 cell = Cell.ascell(cell) 

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

539 special_points=special_points, eps=eps) 

540 

541 

542DEFAULT_KPTS_DENSITY = 5 # points per 1/Angstrom 

543 

544 

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

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

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

548 points = np.concatenate(paths) 

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

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

551 

552 i = 0 

553 for path in paths[:-1]: 

554 i += len(path) 

555 lengths[i - 1] = 0 

556 

557 length = sum(lengths) 

558 

559 if npoints is None: 

560 if density is None: 

561 density = DEFAULT_KPTS_DENSITY 

562 # Set npoints using the length of the path 

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

564 

565 kpts = [] 

566 x0 = 0 

567 x = [] 

568 X = [0] 

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

570 diff = length - x0 

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

572 n = 0 

573 else: 

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

575 

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

577 kpts.append(P + t * d) 

578 x.append(x0 + t * L) 

579 x0 += L 

580 X.append(x0) 

581 if len(points): 

582 kpts.append(points[-1]) 

583 x.append(x0) 

584 

585 if len(kpts) == 0: 

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

587 

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

589 

590 

591get_bandpath = bandpath # old name 

592 

593 

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

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

596 # XXX Should use the Cartesian kpoints. 

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

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

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

600 N = len(kpts) 

601 indices = [] 

602 if N > 0: 

603 indices.append(0) 

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

605 indices.append(N - 1) 

606 return indices 

607 

608 

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

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

611 

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

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

614 as xticklabels. 

615 

616 Parameters: 

617 

618 kpts: list 

619 List of scaled k-points. 

620 

621 cell: list 

622 Unit cell of the atomic structure. 

623 

624 Returns: 

625 

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

627 the second is x coordinates of the special points, 

628 the third is the special points as strings. 

629 """ 

630 

631 if special_points is None: 

632 special_points = get_special_points(cell) 

633 points = np.asarray(kpts) 

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

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

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

637 

638 labels = [] 

639 for kpt in points[indices]: 

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

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

642 break 

643 else: 

644 # No exact match. Try modulus 1: 

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

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

647 break 

648 else: 

649 label = '?' 

650 labels.append(label) 

651 

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

653 return xcoords, ixcoords, labels 

654 

655 

656def indices_to_axis_coords(indices, points, cell): 

657 jump = False # marks a discontinuity in the path 

658 xcoords = [0] 

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

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

661 length = 0 

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

663 else: 

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

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

666 jump = False 

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

668 

669 xcoords = np.array(xcoords) 

670 return xcoords, xcoords[indices] 

671 

672 

673special_paths = { 

674 'cubic': 'GXMGRX,MR', 

675 'fcc': 'GXWKGLUWLK,UX', 

676 'bcc': 'GHNGPH,PN', 

677 'tetragonal': 'GXMGZRAZXR,MA', 

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

679 'hexagonal': 'GMKGALHA,LM,KH', 

680 'monoclinic': 'GYHCEM1AXH1,MDZ,YD', 

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

682 'rhombohedral type 2': 'GPZQGFP1Q1LZ'} 

683 

684 

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

686 """Return dict of special points. 

687 

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

689 Curtarolo:: 

690 

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

692 

693 cell: 3x3 ndarray 

694 Unit cell. 

695 lattice: str 

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

697 bcc, orthorhombic, tetragonal, hexagonal or monoclinic. 

698 eps: float 

699 Tolerance for cell-check. 

700 """ 

701 

702 if isinstance(cell, str): 

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

704 'argument') 

705 lattice, cell = cell, lattice 

706 

707 cell = Cell.ascell(cell) 

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

709 # from the canonical cell to the given one. 

710 # 

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

712 path = cell.bandpath(npoints=0) 

713 return path.special_points 

714 

715 

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

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

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

719 

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

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

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

723 in `path`. 

724 

725 Note 

726 ---- 

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

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

729 

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

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

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

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

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

735 

736 Parameters 

737 ---------- 

738 path: (nk, 3) array-like 

739 Desired path in units of reciprocal lattice vectors. 

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

741 Values on Monkhorst-Pack grid. 

742 icell: (3, 3) array-like 

743 Reciprocal lattice vectors. 

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

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

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

747 Size of Monkhorst-Pack grid. 

748 offset: (3,) array-like 

749 Offset of Monkhorst-Pack grid. 

750 pad_width: int 

751 Padding width to aid interpolation 

752 

753 Returns 

754 ------- 

755 (nbz,) array-like 

756 *values* interpolated to *path*. 

757 """ 

758 from scipy.interpolate import LinearNDInterpolator 

759 

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

761 path = np.dot(path, icell) 

762 

763 # Fold out values from IBZ to BZ: 

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

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

766 

767 # Create padded Monkhorst-Pack grid: 

768 size = np.asarray(size) 

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

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

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

772 k = np.dot(k, icell) 

773 

774 # Fill in boundary values: 

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

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

777 

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

779 interpolated_points = interpolate(path) 

780 

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

782 # try increasing padding 

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

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

785 

786 return interpolated_points 

787 

788 

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

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

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

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

793 

794cc6_1x1 = np.array([ 

795 1, 1, 0, 1, 0, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0, 

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

797 

798cc12_2x3 = np.array([ 

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

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

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

802 

803cc18_sq3xsq3 = np.array([ 

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

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

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

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

808 

809cc18_1x1 = np.array([ 

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

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

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

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

814 

815cc54_sq3xsq3 = np.array([ 

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

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

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