r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import re

2import warnings

3from typing import Dict

5import numpy as np

7import ase # Annotations

8from ase.utils import jsonable

9from ase.cell import Cell

12def monkhorst_pack(size):

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

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

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

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

17 return (kpts + 0.5) / size - 0.5

20def get_monkhorst_pack_size_and_offset(kpts):

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

23 Returns (size, offset), where::

25 kpts = monkhorst_pack(size) + offset.

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

29 if len(kpts) == 1:

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

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])))

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

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

44 kpts0 = monkhorst_pack(size)

45 offsets = kpts - kpts0

47 # All offsets must be identical:

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

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

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

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.

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.

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

71 """

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

73 min_distance, maxperdim, even)

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

77 from ase import Atoms

78 from ase.neighborlist import NeighborList

80 step = 2 if even else 1

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

82 self_interaction=False, bothways=False)

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

88 def generate_mpgrids():

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

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

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

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

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.')

101def get_monkhorst_shape(kpts):

103 return get_monkhorst_pack_size_and_offset(kpts)[0]

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

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

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

110 coordinates, the other is determined.

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.')

124def parse_path_string(s):

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

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.

131 Examples:

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

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

154def resolve_custom_points(pathspec, special_points, eps):

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

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

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

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.

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

169 special_points = dict(special_points)

171 if len(pathspec) == 0:

172 return '', special_points

174 if isinstance(pathspec, str):

175 pathspec = parse_path_string(pathspec)

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,)

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

197 def name_generator():

198 counter = 0

199 while True:

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

201 yield name

202 counter += 1

203 custom_names = name_generator()

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

216 kpt = np.asarray(kpt, float)

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

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

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

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

222 labelseq.append(key)

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(',')

233 last = labelseq.pop()

234 assert last == ','

235 return ''.join(labelseq), special_points

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

249@jsonable('bandpath')

250class BandPath:

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

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:

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])

263 Band paths support JSON I/O:

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

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

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

270 """

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

272 special_points=None, path=None):

273 if kpts is None:

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

276 if special_points is None:

277 special_points = {}

278 else:

279 special_points = normalize_special_points(special_points)

281 if path is None:

282 path = ''

284 cell = Cell(cell)

285 self._cell = cell

286 kpts = np.asarray(kpts)

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

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

289 self._kpts = kpts

290 self._special_points = special_points

291 if not isinstance(path, str):

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

293 self._path = path

295 @property

296 def cell(self) -> Cell:

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

298 return self._cell

300 @property

301 def icell(self) -> Cell:

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

303 return self._icell

305 @property

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

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

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

310 return self._kpts

312 @property

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

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

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

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

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

320 consider deepcopying it."""

321 return self._special_points

323 @property

324 def path(self) -> str:

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

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

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

330 return self._path

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

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

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

336 The operation will typically be a permutation/flipping

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

338 # XXX acceptable operations are probably only those

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

340 #

341 # We should insert a check.

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

343 # if they change lengths, volume etc.

344 special_points = {}

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

346 special_points[name] = value @ op

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

349 special_points=special_points,

350 path=self.path)

352 def todict(self):

353 return {'kpts': self.kpts,

354 'special_points': self.special_points,

355 'labelseq': self.path,

356 'cell': self.cell}

358 def interpolate(

359 self,

360 path: str = None,

361 npoints: int = None,

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

363 density: float = None,

364 ) -> 'BandPath':

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

366 if path is None:

367 path = self.path

369 if special_points is None:

370 special_points = self.special_points

372 pathnames, pathcoords = resolve_kpt_path_string(path, special_points)

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

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

375 special_points=special_points)

377 def _scale(self, coords):

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

380 def __repr__(self):

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

382 .format(self.__class__.__name__,

383 repr(self.path),

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

385 len(self.kpts)))

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

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

389 return self._scale(self.kpts)

391 def __iter__(self):

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

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

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

396 special_x_coords), and people would use tuple unpacking for

397 those.

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

400 It will be removed in the future.

402 """

403 import warnings

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

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

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

407 yield self.kpts

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

410 yield x

411 yield xspecial

413 def __getitem__(self, index):

414 # Temp compatibility stuff, see __iter__

415 return tuple(self)[index]

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

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

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

422 index2name = self._find_special_point_indices(eps)

423 indices = sorted(index2name)

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

425 xcoords, special_xcoords = indices_to_axis_coords(

426 indices, self.kpts, self.cell)

427 return xcoords, special_xcoords, labels

429 def _find_special_point_indices(self, eps):

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

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

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

434 # to fully make sense!

435 index2name = {}

436 nkpts = len(self.kpts)

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

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

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

441 args = np.argwhere(distances < eps)

442 for arg in args.flat:

443 dist = distances[arg]

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

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

446 if not any(neighbours < dist):

447 index2name[arg] = name

449 return index2name

451 def plot(self, **plotkwargs):

452 """Visualize this bandpath.

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

455 import ase.dft.bz as bz

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

458 plotkwargs.pop('dimension', None)

460 special_points = self.special_points

461 labelseq, coords = resolve_kpt_path_string(self.path,

462 special_points)

464 paths = []

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

467 subpath_coords = np.array(subpath_coords)

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

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

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

474 if label not in points_already_plotted:

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

477 kw = {'vectors': True,

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

480 kw.update(plotkwargs)

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

482 points=self.cartesian_kpts(),

483 **kw)

485 def free_electron_band_structure(

486 self, **kwargs

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

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

490 Keyword arguments are passed to

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

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

494 from ase import Atoms

495 from ase.calculators.test import FreeElectrons

496 from ase.spectrum.band_structure import calculate_band_structure

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

498 atoms.calc = FreeElectrons(**kwargs)

499 bs = calculate_band_structure(atoms, path=self)

500 return bs

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

504 eps=2e-4):

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

507 path: list or str

508 Can be:

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

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

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

513 cell: 3x3

514 Unit cell of the atoms.

515 npoints: int

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

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

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

519 density: float

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

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

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

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

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

525 will be used.

526 special_points: dict or None

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

528 points will be derived from the cell.

529 eps: float

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

532 You may define npoints or density but not both.

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

536 cell = Cell.ascell(cell)

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

538 special_points=special_points, eps=eps)

541DEFAULT_KPTS_DENSITY = 5 # points per 1/Angstrom

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

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

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

547 points = np.concatenate(paths)

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

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

551 i = 0

552 for path in paths[:-1]:

553 i += len(path)

554 lengths[i - 1] = 0

556 length = sum(lengths)

558 if npoints is None:

559 if density is None:

560 density = DEFAULT_KPTS_DENSITY

561 # Set npoints using the length of the path

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

564 kpts = []

565 x0 = 0

566 x = []

567 X = [0]

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

569 diff = length - x0

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

571 n = 0

572 else:

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

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

576 kpts.append(P + t * d)

577 x.append(x0 + t * L)

578 x0 += L

579 X.append(x0)

580 if len(points):

581 kpts.append(points[-1])

582 x.append(x0)

584 if len(kpts) == 0:

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

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

590get_bandpath = bandpath # old name

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

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

595 # XXX Should use the Cartesian kpoints.

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

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

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

599 N = len(kpts)

600 indices = []

601 if N > 0:

602 indices.append(0)

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

604 indices.append(N - 1)

605 return indices

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

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

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

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

613 as xticklabels.

615 Parameters:

617 kpts: list

618 List of scaled k-points.

620 cell: list

621 Unit cell of the atomic structure.

623 Returns:

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

626 the second is x coordinates of the special points,

627 the third is the special points as strings.

628 """

630 if special_points is None:

631 special_points = get_special_points(cell)

632 points = np.asarray(kpts)

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

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

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

637 labels = []

638 for kpt in points[indices]:

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

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

641 break

642 else:

643 # No exact match. Try modulus 1:

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

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

646 break

647 else:

648 label = '?'

649 labels.append(label)

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

652 return xcoords, ixcoords, labels

655def indices_to_axis_coords(indices, points, cell):

656 jump = False # marks a discontinuity in the path

657 xcoords = [0]

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

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

660 length = 0

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

662 else:

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

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

665 jump = False

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

668 xcoords = np.array(xcoords)

669 return xcoords, xcoords[indices]

672special_paths = {

673 'cubic': 'GXMGRX,MR',

674 'fcc': 'GXWKGLUWLK,UX',

675 'bcc': 'GHNGPH,PN',

676 'tetragonal': 'GXMGZRAZXR,MA',

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

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

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

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

681 'rhombohedral type 2': 'GPZQGFP1Q1LZ'}

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

685 """Return dict of special points.

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

688 Curtarolo::

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

692 cell: 3x3 ndarray

693 Unit cell.

694 lattice: str

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

696 bcc, orthorhombic, tetragonal, hexagonal or monoclinic.

697 eps: float

698 Tolerance for cell-check.

699 """

701 if isinstance(cell, str):

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

703 'argument')

704 lattice, cell = cell, lattice

706 cell = Cell.ascell(cell)

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

708 # from the canonical cell to the given one.

709 #

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

711 path = cell.bandpath(npoints=0)

712 return path.special_points

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

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

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

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

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

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

722 in `path`.

724 Note

725 ----

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

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

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

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

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

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

735 Parameters

736 ----------

737 path: (nk, 3) array-like

738 Desired path in units of reciprocal lattice vectors.

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

740 Values on Monkhorst-Pack grid.

741 icell: (3, 3) array-like

742 Reciprocal lattice vectors.

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

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

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

746 Size of Monkhorst-Pack grid.

747 offset: (3,) array-like

748 Offset of Monkhorst-Pack grid.

750 Padding width to aid interpolation

752 Returns

753 -------

754 (nbz,) array-like

755 *values* interpolated to *path*.

756 """

757 from scipy.interpolate import LinearNDInterpolator

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

760 path = np.dot(path, icell)

762 # Fold out values from IBZ to BZ:

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

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

766 # Create padded Monkhorst-Pack grid:

767 size = np.asarray(size)

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

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

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

771 k = np.dot(k, icell)

773 # Fill in boundary values:

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

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

778 interpolated_points = interpolate(path)

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

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

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

785 return interpolated_points

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

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

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

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

793cc6_1x1 = np.array([

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

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

797cc12_2x3 = np.array([

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

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

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

802cc18_sq3xsq3 = np.array([

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

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

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

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

808cc18_1x1 = np.array([

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

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

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

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

814cc54_sq3xsq3 = np.array([

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

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

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

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

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

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

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

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

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

825cc54_1x1 = np.array([

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

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

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

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

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

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

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

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

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

836cc162_sq3xsq3 = np.array([

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

865cc162_1x1 = np.array([

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

898#

899# In units of the reciprocal basis vectors.

900#

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

902sc_special_points = {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

939# Old version of dictionary kept for backwards compatibility.

940# Not for ordinary use.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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