Coverage for /builds/ase/ase/ase/dft/kpoints.py: 86.59%

1# fmt: off 

2 

3import re 

4import warnings 

5from typing import Dict 

6 

7import numpy as np 

8 

9import ase # Annotations 

10from ase.cell import Cell 

11from ase.utils import jsonable 

12 

13 

14def monkhorst_pack(size): 

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

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

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

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

19 return (kpts + 0.5) / size - 0.5 

20 

21 

22def get_monkhorst_pack_size_and_offset(kpts): 

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

24 

25 Returns (size, offset), where:: 

26 

27 kpts = monkhorst_pack(size) + offset. 

28 

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

30 

31 if len(kpts) == 1: 

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

33 

34 size = np.zeros(3, int) 

35 for c in range(3): 

36 # Determine increment between k-points along current axis 

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

38 

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

40 if delta > 1e-8: 

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

42 else: 

43 size[c] = 1 

44 

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

46 kpts0 = monkhorst_pack(size) 

47 offsets = kpts - kpts0 

48 

49 # All offsets must be identical: 

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

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

52 

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

54 

55 

56def mindistance2monkhorstpack(atoms, *, 

57 min_distance, 

58 maxperdim=16, 

59 even=True): 

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

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

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

63 (nx, ny, nz)-supercell. 

64 

65 Compared to ase.calculators.calculator kptdensity2monkhorstpack 

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

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

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

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

70 invariant to particular choice of cell representations. 

71 

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

73 """ 

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

75 min_distance, maxperdim, even) 

76 

77 

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

79 from ase import Atoms 

80 from ase.neighborlist import NeighborList 

81 

82 step = 2 if even else 1 

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

84 self_interaction=False, bothways=False) 

85 

86 def check(nkpts_c): 

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

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

89 

90 def generate_mpgrids(): 

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

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

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

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

95 

96 try: 

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

98 except StopIteration: 

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

100 ' Try running with a larger maxperdim.') 

101 

102 

103def get_monkhorst_shape(kpts): 

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

105 return get_monkhorst_pack_size_and_offset(kpts)[0] 

106 

107 

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

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

110 

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

112 coordinates, the other is determined. 

113 

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

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

116 """ 

117 if ckpts_kv is None: 

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

119 return np.dot(skpts_kc, icell_cv) 

120 elif skpts_kc is None: 

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

122 else: 

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

124 

125 

126def parse_path_string(s): 

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

128 

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

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

131 list of labels. 

132 

133 Examples: 

134 

135 >>> parse_path_string('GX') 

136 [['G', 'X']] 

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

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

139 """ 

140 paths = [] 

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

142 names = [name 

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

144 if name] 

145 paths.append(names) 

146 return paths 

147 

148 

149def resolve_kpt_path_string(path, special_points): 

150 paths = parse_path_string(path) 

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

152 for subpath in paths] 

153 return paths, coords 

154 

155 

156def resolve_custom_points(pathspec, special_points, eps): 

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

158 

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

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

161 

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

163 are generated dynamically as necessary, updating the special_point 

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

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

166 labelled as such.""" 

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

168 # be lazy and call it on scaled ones. 

169 

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

171 special_points = dict(special_points) 

172 

173 if len(pathspec) == 0: 

174 return '', special_points 

175 

176 if isinstance(pathspec, str): 

177 pathspec = parse_path_string(pathspec) 

178 

179 def looks_like_single_kpoint(obj): 

180 if isinstance(obj, str): 

181 return True 

182 try: 

183 arr = np.asarray(obj, float) 

184 except ValueError: 

185 return False 

186 else: 

187 return arr.shape == (3,) 

188 

189 # We accept inputs that are either 

190 # 1) string notation 

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

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

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

194 for thing in pathspec: 

195 if looks_like_single_kpoint(thing): 

196 pathspec = [pathspec] 

197 break 

198 

199 def name_generator(): 

200 counter = 0 

201 while True: 

202 name = f'Kpt{counter}' 

203 yield name 

204 counter += 1 

205 

206 custom_names = name_generator() 

207 

208 labelseq = [] 

209 for subpath in pathspec: 

210 for kpt in subpath: 

211 if isinstance(kpt, str): 

212 if kpt not in special_points: 

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

214 .format(kpt, 

215 ''.join(special_points))) 

216 labelseq.append(kpt) 

217 continue 

218 

219 kpt = np.asarray(kpt, float) 

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

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

222 

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

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

225 labelseq.append(key) 

226 break # Already present 

227 else: 

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

229 name = next(custom_names) 

230 while name in special_points: 

231 name = next(custom_names) 

232 special_points[name] = kpt 

233 labelseq.append(name) 

234 labelseq.append(',') 

235 

236 last = labelseq.pop() 

237 assert last == ',' 

238 return ''.join(labelseq), special_points 

239 

240 

241def normalize_special_points(special_points): 

242 dct = {} 

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

244 if not isinstance(name, str): 

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

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

247 raise ValueError('Expected 3 kpoint coordinates') 

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

249 return dct 

250 

251 

252@jsonable('bandpath') 

253class BandPath: 

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

255 

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

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

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

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

260 

261 >>> from ase.lattice import CUB 

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

263 >>> path 

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

265 

266 Band paths support JSON I/O: 

267 

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

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

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

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

272 

273 """ 

274 

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

276 special_points=None, path=None): 

277 if kpts is None: 

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

279 

280 if special_points is None: 

281 special_points = {} 

282 else: 

283 special_points = normalize_special_points(special_points) 

284 

285 if path is None: 

286 path = '' 

287 

288 cell = Cell(cell) 

289 self._cell = cell 

290 kpts = np.asarray(kpts) 

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

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

293 self._kpts = kpts 

294 self._special_points = special_points 

295 if not isinstance(path, str): 

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

297 self._path = path 

298 

299 @property 

300 def cell(self) -> Cell: 

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

302 return self._cell 

303 

304 @property 

305 def icell(self) -> Cell: 

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

307 return self._icell 

308 

309 @property 

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

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

312 

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

314 return self._kpts 

315 

316 @property 

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

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

319 

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

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

322 

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

324 consider deepcopying it.""" 

325 return self._special_points 

326 

327 @property 

328 def path(self) -> str: 

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

330 

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

332 

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

334 return self._path 

335 

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

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

338 

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

340 The operation will typically be a permutation/flipping 

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

342 # XXX acceptable operations are probably only those 

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

344 # 

345 # We should insert a check. 

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

347 # if they change lengths, volume etc. 

348 special_points = { 

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

350 } 

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

352 special_points=special_points, 

353 path=self.path) 

354 

355 def todict(self): 

356 return {'kpts': self.kpts, 

357 'special_points': self.special_points, 

358 'labelseq': self.path, 

359 'cell': self.cell} 

360 

361 def interpolate( 

362 self, 

363 path: str = None, 

364 npoints: int = None, 

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

366 density: float = None, 

367 ) -> 'BandPath': 

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

369 if path is None: 

370 path = self.path 

371 

372 if special_points is None: 

373 special_points = self.special_points 

374 

375 _pathnames, pathcoords = resolve_kpt_path_string(path, special_points) 

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

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

378 special_points=special_points) 

379 

380 def _scale(self, coords): 

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

382 

383 def __repr__(self): 

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

385 .format(self.__class__.__name__, 

386 repr(self.path), 

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

388 len(self.kpts))) 

389 

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

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

392 return self._scale(self.kpts) 

393 

394 def __iter__(self): 

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

396 

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

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

399 special_x_coords), and people would use tuple unpacking for 

400 those. 

401 

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

403 It will be removed in the future. 

404 

405 """ 

406 import warnings 

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

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

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

410 yield self.kpts 

411 

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

413 yield x 

414 yield xspecial 

415 

416 def __getitem__(self, index): 

417 # Temp compatibility stuff, see __iter__ 

418 return tuple(self)[index] 

419 

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

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

422 

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

424 

425 index2name = self._find_special_point_indices(eps) 

426 indices = sorted(index2name) 

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

428 xcoords, special_xcoords = indices_to_axis_coords( 

429 indices, self.kpts, self.cell) 

430 return xcoords, special_xcoords, labels 

431 

432 def _find_special_point_indices(self, eps): 

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

434 

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

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

437 # to fully make sense! 

438 index2name = {} 

439 nkpts = len(self.kpts) 

440 

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

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

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

444 args = np.argwhere(distances < eps) 

445 for arg in args.flat: 

446 dist = distances[arg] 

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

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

449 if not any(neighbours < dist): 

450 index2name[arg] = name 

451 

452 return index2name 

453 

454 def plot(self, **plotkwargs): 

455 """Visualize this bandpath. 

456 

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

458 import ase.dft.bz as bz 

459 

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

461 plotkwargs.pop('dimension', None) 

462 

463 special_points = self.special_points 

464 labelseq, coords = resolve_kpt_path_string(self.path, 

465 special_points) 

466 

467 paths = [] 

468 points_already_plotted = set() 

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

470 subpath_coords = np.array(subpath_coords) 

471 points_already_plotted.update(subpath_labels) 

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

473 

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

475 # not plotted already: 

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

477 if label not in points_already_plotted: 

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

479 

480 kw = {'vectors': True, 

481 'pointstyle': {'marker': '.'}} 

482 

483 kw.update(plotkwargs) 

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

485 points=self.cartesian_kpts(), 

486 **kw) 

487 

488 def free_electron_band_structure( 

489 self, **kwargs 

490 ) -> 'ase.spectrum.band_structure.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 import Atoms 

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 kpts: list 

621 List of scaled k-points. 

622 

623 cell: list 

624 Unit cell of the atomic structure. 

625 

626 Returns: 

627 

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

629 the second is x coordinates of the special points, 

630 the third is the special points as strings. 

631 """ 

632 

633 if special_points is None: 

634 special_points = get_special_points(cell) 

635 points = np.asarray(kpts) 

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

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

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

639 

640 labels = [] 

641 for kpt in points[indices]: 

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

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

644 break 

645 else: 

646 # No exact match. Try modulus 1: 

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

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

649 break 

650 else: 

651 label = '?' 

652 labels.append(label) 

653 

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

655 return xcoords, ixcoords, labels 

656 

657 

658def indices_to_axis_coords(indices, points, cell): 

659 jump = False # marks a discontinuity in the path 

660 xcoords = [0] 

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

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

663 length = 0 

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

665 else: 

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

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

668 jump = False 

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

670 

671 xcoords = np.array(xcoords) 

672 return xcoords, xcoords[indices] 

673 

674 

675special_paths = { 

676 'cubic': 'GXMGRX,MR', 

677 'fcc': 'GXWKGLUWLK,UX', 

678 'bcc': 'GHNGPH,PN', 

679 'tetragonal': 'GXMGZRAZXR,MA', 

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

681 'hexagonal': 'GMKGALHA,LM,KH', 

682 'monoclinic': 'GYHCEM1AXH1,MDZ,YD', 

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

684 'rhombohedral type 2': 'GPZQGFP1Q1LZ'} 

685 

686 

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

688 """Return dict of special points. 

689 

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

691 Curtarolo:: 

692 

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

694 

695 cell: 3x3 ndarray 

696 Unit cell. 

697 lattice: str 

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

699 bcc, orthorhombic, tetragonal, hexagonal or monoclinic. 

700 eps: float 

701 Tolerance for cell-check. 

702 """ 

703 

704 if isinstance(cell, str): 

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

706 'argument') 

707 lattice, cell = cell, lattice 

708 

709 cell = Cell.ascell(cell) 

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

711 # from the canonical cell to the given one. 

712 # 

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

714 path = cell.bandpath(npoints=0) 

715 return path.special_points 

716 

717 

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

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

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

721 

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

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

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

725 in `path`. 

726 

727 Note 

728 ---- 

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

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

731 

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

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

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

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

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

737 

738 Parameters 

739 ---------- 

740 path: (nk, 3) array-like 

741 Desired path in units of reciprocal lattice vectors. 

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

743 Values on Monkhorst-Pack grid. 

744 icell: (3, 3) array-like 

745 Reciprocal lattice vectors. 

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

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

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

749 Size of Monkhorst-Pack grid. 

750 offset: (3,) array-like 

751 Offset of Monkhorst-Pack grid. 

752 pad_width: int 

753 Padding width to aid interpolation 

754 

755 Returns 

756 ------- 

757 (nbz,) array-like 

758 *values* interpolated to *path*. 

759 """ 

760 from scipy.interpolate import LinearNDInterpolator 

761 

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

763 path = np.dot(path, icell) 

764 

765 # Fold out values from IBZ to BZ: 

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

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

768 

769 # Create padded Monkhorst-Pack grid: 

770 size = np.asarray(size) 

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

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

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

774 k = np.dot(k, icell) 

775 

776 # Fill in boundary values: 

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

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

779 

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

781 interpolated_points = interpolate(path) 

782 

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

784 # try increasing padding 

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

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

787 

788 return interpolated_points 

789 

790 

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

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

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

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

795 

796cc6_1x1 = np.array([ 

797 1, 1, 0, 1, 0, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0, 

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

799 

800cc12_2x3 = np.array([ 

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

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

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

804 

805cc18_sq3xsq3 = np.array([ 

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

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

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

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

810 

811cc18_1x1 = np.array([ 

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

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

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

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

816 

817cc54_sq3xsq3 = np.array([ 

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

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

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