Coverage for /builds/debichem-team/python-ase/ase/dft/kpoints.py : 83.44%
Hot-keys on this page
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 get_monkhorst_shape(kpts):
55 warnings.warn('Use get_monkhorst_pack_size_and_offset()[0] instead.')
56 return get_monkhorst_pack_size_and_offset(kpts)[0]
59def kpoint_convert(cell_cv, skpts_kc=None, ckpts_kv=None):
60 """Convert k-points between scaled and cartesian coordinates.
62 Given the atomic unit cell, and either the scaled or cartesian k-point
63 coordinates, the other is determined.
65 The k-point arrays can be either a single point, or a list of points,
66 i.e. the dimension k can be empty or multidimensional.
67 """
68 if ckpts_kv is None:
69 icell_cv = 2 * np.pi * np.linalg.pinv(cell_cv).T
70 return np.dot(skpts_kc, icell_cv)
71 elif skpts_kc is None:
72 return np.dot(ckpts_kv, cell_cv.T) / (2 * np.pi)
73 else:
74 raise KeyError('Either scaled or cartesian coordinates must be given.')
77def parse_path_string(s):
78 """Parse compact string representation of BZ path.
80 A path string can have several non-connected sections separated by
81 commas. The return value is a list of sections where each section is a
82 list of labels.
84 Examples:
86 >>> parse_path_string('GX')
87 [['G', 'X']]
88 >>> parse_path_string('GX,M1A')
89 [['G', 'X'], ['M1', 'A']]
90 """
91 paths = []
92 for path in s.split(','):
93 names = [name
94 for name in re.split(r'([A-Z][a-z0-9]*)', path)
95 if name]
96 paths.append(names)
97 return paths
100def resolve_kpt_path_string(path, special_points):
101 paths = parse_path_string(path)
102 coords = [np.array([special_points[sym] for sym in subpath]).reshape(-1, 3)
103 for subpath in paths]
104 return paths, coords
107def resolve_custom_points(pathspec, special_points, eps):
108 """Resolve a path specification into a string.
110 The path specification is a list path segments, each segment being a kpoint
111 label or kpoint coordinate, or a single such segment.
113 Return a string representing the same path. Generic kpoint labels
114 are generated dynamically as necessary, updating the special_point
115 dictionary if necessary. The tolerance eps is used to see whether
116 coordinates are close enough to a special point to deserve being
117 labelled as such."""
118 # This should really run on Cartesian coordinates but we'll probably
119 # be lazy and call it on scaled ones.
121 # We may add new points below so take a copy of the input:
122 special_points = dict(special_points)
124 if len(pathspec) == 0:
125 return '', special_points
127 if isinstance(pathspec, str):
128 pathspec = parse_path_string(pathspec)
130 def looks_like_single_kpoint(obj):
131 if isinstance(obj, str):
132 return True
133 try:
134 arr = np.asarray(obj, float)
135 except ValueError:
136 return False
137 else:
138 return arr.shape == (3,)
140 # We accept inputs that are either
141 # 1) string notation
142 # 2) list of kpoints (each either a string or three floats)
143 # 3) list of list of kpoints; each toplevel list is a contiguous subpath
144 # Here we detect form 2 and normalize to form 3:
145 for thing in pathspec:
146 if looks_like_single_kpoint(thing):
147 pathspec = [pathspec]
148 break
150 def name_generator():
151 counter = 0
152 while True:
153 name = 'Kpt{}'.format(counter)
154 yield name
155 counter += 1
156 custom_names = name_generator()
158 labelseq = []
159 for subpath in pathspec:
160 for kpt in subpath:
161 if isinstance(kpt, str):
162 if kpt not in special_points:
163 raise KeyError('No kpoint "{}" among "{}"'
164 .format(kpt,
165 ''.join(special_points)))
166 labelseq.append(kpt)
167 continue
169 kpt = np.asarray(kpt, float)
170 if not kpt.shape == (3,):
171 raise ValueError(f'Not a valid kpoint: {kpt}')
173 for key, val in special_points.items():
174 if np.abs(kpt - val).max() < eps:
175 labelseq.append(key)
176 break # Already present
177 else:
178 # New special point - search for name we haven't used yet:
179 name = next(custom_names)
180 while name in special_points:
181 name = next(custom_names)
182 special_points[name] = kpt
183 labelseq.append(name)
184 labelseq.append(',')
186 last = labelseq.pop()
187 assert last == ','
188 return ''.join(labelseq), special_points
191def normalize_special_points(special_points):
192 dct = {}
193 for name, value in special_points.items():
194 if not isinstance(name, str):
195 raise TypeError('Expected name to be a string')
196 if not np.shape(value) == (3,):
197 raise ValueError('Expected 3 kpoint coordinates')
198 dct[name] = np.asarray(value, float)
199 return dct
202@jsonable('bandpath')
203class BandPath:
204 """Represents a Brillouin zone path or bandpath.
206 A band path has a unit cell, a path specification, special points,
207 and interpolated k-points. Band paths are typically created
208 indirectly using the :class:`~ase.geometry.Cell` or
209 :class:`~ase.lattice.BravaisLattice` classes:
211 >>> from ase.lattice import CUB
212 >>> path = CUB(3).bandpath()
213 >>> path
214 BandPath(path='GXMGRX,MR', cell=[3x3], special_points={GMRX}, kpts=[40x3])
216 Band paths support JSON I/O:
218 >>> from ase.io.jsonio import read_json
219 >>> path.write('mybandpath.json')
220 >>> read_json('mybandpath.json')
221 BandPath(path='GXMGRX,MR', cell=[3x3], special_points={GMRX}, kpts=[40x3])
223 """
224 def __init__(self, cell, kpts=None,
225 special_points=None, path=None):
226 if kpts is None:
227 kpts = np.empty((0, 3))
229 if special_points is None:
230 special_points = {}
231 else:
232 special_points = normalize_special_points(special_points)
234 if path is None:
235 path = ''
237 cell = Cell(cell)
238 self._cell = cell
239 kpts = np.asarray(kpts)
240 assert kpts.ndim == 2 and kpts.shape[1] == 3 and kpts.dtype == float
241 self._icell = self.cell.reciprocal()
242 self._kpts = kpts
243 self._special_points = special_points
244 if not isinstance(path, str):
245 raise TypeError(f'path must be a string; was {path!r}')
246 self._path = path
248 @property
249 def cell(self) -> Cell:
250 """The :class:`~ase.cell.Cell` of this BandPath."""
251 return self._cell
253 @property
254 def icell(self) -> Cell:
255 """Reciprocal cell of this BandPath as a :class:`~ase.cell.Cell`."""
256 return self._icell
258 @property
259 def kpts(self) -> np.ndarray:
260 """The kpoints of this BandPath as an array of shape (nkpts, 3).
262 The kpoints are given in units of the reciprocal cell."""
263 return self._kpts
265 @property
266 def special_points(self) -> Dict[str, np.ndarray]:
267 """Special points of this BandPath as a dictionary.
269 The dictionary maps names (such as `'G'`) to kpoint coordinates
270 in units of the reciprocal cell as a 3-element numpy array.
272 It's unwise to edit this dictionary directly. If you need that,
273 consider deepcopying it."""
274 return self._special_points
276 @property
277 def path(self) -> str:
278 """The string specification of this band path.
280 This is a specification of the form `'GXWKGLUWLK,UX'`.
282 Comma marks a discontinuous jump: K is not connected to U."""
283 return self._path
285 def transform(self, op: np.ndarray) -> 'BandPath':
286 """Apply 3x3 matrix to this BandPath and return new BandPath.
288 This is useful for converting the band path to another cell.
289 The operation will typically be a permutation/flipping
290 established by a function such as Niggli reduction."""
291 # XXX acceptable operations are probably only those
292 # who come from Niggli reductions (permutations etc.)
293 #
294 # We should insert a check.
295 # I wonder which operations are valid? They won't be valid
296 # if they change lengths, volume etc.
297 special_points = {}
298 for name, value in self.special_points.items():
299 special_points[name] = value @ op
301 return BandPath(op.T @ self.cell, kpts=self.kpts @ op,
302 special_points=special_points,
303 path=self.path)
305 def todict(self):
306 return {'kpts': self.kpts,
307 'special_points': self.special_points,
308 'labelseq': self.path,
309 'cell': self.cell}
311 def interpolate(
312 self,
313 path: str = None,
314 npoints: int = None,
315 special_points: Dict[str, np.ndarray] = None,
316 density: float = None,
317 ) -> 'BandPath':
318 """Create new bandpath, (re-)interpolating kpoints from this one."""
319 if path is None:
320 path = self.path
322 if special_points is None:
323 special_points = self.special_points
325 pathnames, pathcoords = resolve_kpt_path_string(path, special_points)
326 kpts, x, X = paths2kpts(pathcoords, self.cell, npoints, density)
327 return BandPath(self.cell, kpts, path=path,
328 special_points=special_points)
330 def _scale(self, coords):
331 return np.dot(coords, self.icell)
333 def __repr__(self):
334 return ('{}(path={}, cell=[3x3], special_points={{{}}}, kpts=[{}x3])'
335 .format(self.__class__.__name__,
336 repr(self.path),
337 ''.join(sorted(self.special_points)),
338 len(self.kpts)))
340 def cartesian_kpts(self) -> np.ndarray:
341 """Get Cartesian kpoints from this bandpath."""
342 return self._scale(self.kpts)
344 def __iter__(self):
345 """XXX Compatibility hack for bandpath() function.
347 bandpath() now returns a BandPath object, which is a Good
348 Thing. However it used to return a tuple of (kpts, x_axis,
349 special_x_coords), and people would use tuple unpacking for
350 those.
352 This function makes tuple unpacking work in the same way.
353 It will be removed in the future.
355 """
356 import warnings
357 warnings.warn('Please do not use (kpts, x, X) = bandpath(...). '
358 'Use path = bandpath(...) and then kpts = path.kpts and '
359 '(x, X, labels) = path.get_linear_kpoint_axis().')
360 yield self.kpts
362 x, xspecial, _ = self.get_linear_kpoint_axis()
363 yield x
364 yield xspecial
366 def __getitem__(self, index):
367 # Temp compatibility stuff, see __iter__
368 return tuple(self)[index]
370 def get_linear_kpoint_axis(self, eps=1e-5):
371 """Define x axis suitable for plotting a band structure.
373 See :func:`ase.dft.kpoints.labels_from_kpts`."""
375 index2name = self._find_special_point_indices(eps)
376 indices = sorted(index2name)
377 labels = [index2name[index] for index in indices]
378 xcoords, special_xcoords = indices_to_axis_coords(
379 indices, self.kpts, self.cell)
380 return xcoords, special_xcoords, labels
382 def _find_special_point_indices(self, eps):
383 """Find indices of kpoints which are close to special points.
385 The result is returned as a dictionary mapping indices to labels."""
386 # XXX must work in Cartesian coordinates for comparison to eps
387 # to fully make sense!
388 index2name = {}
389 nkpts = len(self.kpts)
391 for name, kpt in self.special_points.items():
392 displacements = self.kpts - kpt[np.newaxis, :]
393 distances = np.linalg.norm(displacements, axis=1)
394 args = np.argwhere(distances < eps)
395 for arg in args.flat:
396 dist = distances[arg]
397 # Check if an adjacent point exists and is even closer:
398 neighbours = distances[max(arg - 1, 0):min(arg + 1, nkpts - 1)]
399 if not any(neighbours < dist):
400 index2name[arg] = name
402 return index2name
404 def plot(self, **plotkwargs):
405 """Visualize this bandpath.
407 Plots the irreducible Brillouin zone and this bandpath."""
408 import ase.dft.bz as bz
410 # We previously had a "dimension=3" argument which is now unused.
411 plotkwargs.pop('dimension', None)
413 special_points = self.special_points
414 labelseq, coords = resolve_kpt_path_string(self.path,
415 special_points)
417 paths = []
418 points_already_plotted = set()
419 for subpath_labels, subpath_coords in zip(labelseq, coords):
420 subpath_coords = np.array(subpath_coords)
421 points_already_plotted.update(subpath_labels)
422 paths.append((subpath_labels, self._scale(subpath_coords)))
424 # Add each special point as a single-point subpath if they were
425 # not plotted already:
426 for label, point in special_points.items():
427 if label not in points_already_plotted:
428 paths.append(([label], [self._scale(point)]))
430 kw = {'vectors': True,
431 'pointstyle': {'marker': '.'}}
433 kw.update(plotkwargs)
434 return bz.bz_plot(self.cell, paths=paths,
435 points=self.cartesian_kpts(),
436 **kw)
438 def free_electron_band_structure(
439 self, **kwargs
440 ) -> 'ase.spectrum.band_structure.BandStructure':
441 """Return band structure of free electrons for this bandpath.
443 Keyword arguments are passed to
444 :class:`~ase.calculators.test.FreeElectrons`.
446 This is for mostly testing and visualization."""
447 from ase import Atoms
448 from ase.calculators.test import FreeElectrons
449 from ase.spectrum.band_structure import calculate_band_structure
450 atoms = Atoms(cell=self.cell, pbc=True)
451 atoms.calc = FreeElectrons(**kwargs)
452 bs = calculate_band_structure(atoms, path=self)
453 return bs
456def bandpath(path, cell, npoints=None, density=None, special_points=None,
457 eps=2e-4):
458 """Make a list of kpoints defining the path between the given points.
460 path: list or str
461 Can be:
463 * a string that parse_path_string() understands: 'GXL'
464 * a list of BZ points: [(0, 0, 0), (0.5, 0, 0)]
465 * or several lists of BZ points if the the path is not continuous.
466 cell: 3x3
467 Unit cell of the atoms.
468 npoints: int
469 Length of the output kpts list. If too small, at least the beginning
470 and ending point of each path segment will be used. If None (default),
471 it will be calculated using the supplied density or a default one.
472 density: float
473 k-points per 1/A on the output kpts list. If npoints is None,
474 the number of k-points in the output list will be:
475 npoints = density * path total length (in Angstroms).
476 If density is None (default), use 5 k-points per A⁻¹.
477 If the calculated npoints value is less than 50, a minimum value of 50
478 will be used.
479 special_points: dict or None
480 Dictionary mapping names to special points. If None, the special
481 points will be derived from the cell.
482 eps: float
483 Precision used to identify Bravais lattice, deducing special points.
485 You may define npoints or density but not both.
487 Return a :class:`~ase.dft.kpoints.BandPath` object."""
489 cell = Cell.ascell(cell)
490 return cell.bandpath(path, npoints=npoints, density=density,
491 special_points=special_points, eps=eps)
494DEFAULT_KPTS_DENSITY = 5 # points per 1/Angstrom
497def paths2kpts(paths, cell, npoints=None, density=None):
498 if not(npoints is None or density is None):
499 raise ValueError('You may define npoints or density, but not both.')
500 points = np.concatenate(paths)
501 dists = points[1:] - points[:-1]
502 lengths = [np.linalg.norm(d) for d in kpoint_convert(cell, skpts_kc=dists)]
504 i = 0
505 for path in paths[:-1]:
506 i += len(path)
507 lengths[i - 1] = 0
509 length = sum(lengths)
511 if npoints is None:
512 if density is None:
513 density = DEFAULT_KPTS_DENSITY
514 # Set npoints using the length of the path
515 npoints = int(round(length * density))
517 kpts = []
518 x0 = 0
519 x = []
520 X = [0]
521 for P, d, L in zip(points[:-1], dists, lengths):
522 diff = length - x0
523 if abs(diff) < 1e-6:
524 n = 0
525 else:
526 n = max(2, int(round(L * (npoints - len(x)) / diff)))
528 for t in np.linspace(0, 1, n)[:-1]:
529 kpts.append(P + t * d)
530 x.append(x0 + t * L)
531 x0 += L
532 X.append(x0)
533 if len(points):
534 kpts.append(points[-1])
535 x.append(x0)
537 if len(kpts) == 0:
538 kpts = np.empty((0, 3))
540 return np.array(kpts), np.array(x), np.array(X)
543get_bandpath = bandpath # old name
546def find_bandpath_kinks(cell, kpts, eps=1e-5):
547 """Find indices of those kpoints that are not interior to a line segment."""
548 # XXX Should use the Cartesian kpoints.
549 # Else comparison to eps will be anisotropic and hence arbitrary.
550 diffs = kpts[1:] - kpts[:-1]
551 kinks = abs(diffs[1:] - diffs[:-1]).sum(1) > eps
552 N = len(kpts)
553 indices = []
554 if N > 0:
555 indices.append(0)
556 indices.extend(np.arange(1, N - 1)[kinks])
557 indices.append(N - 1)
558 return indices
561def labels_from_kpts(kpts, cell, eps=1e-5, special_points=None):
562 """Get an x-axis to be used when plotting a band structure.
564 The first of the returned lists can be used as a x-axis when plotting
565 the band structure. The second list can be used as xticks, and the third
566 as xticklabels.
568 Parameters:
570 kpts: list
571 List of scaled k-points.
573 cell: list
574 Unit cell of the atomic structure.
576 Returns:
578 Three arrays; the first is a list of cumulative distances between k-points,
579 the second is x coordinates of the special points,
580 the third is the special points as strings.
581 """
583 if special_points is None:
584 special_points = get_special_points(cell)
585 points = np.asarray(kpts)
586 # XXX Due to this mechanism, we are blind to special points
587 # that lie on straight segments such as [K, G, -K].
588 indices = find_bandpath_kinks(cell, kpts, eps=1e-5)
590 labels = []
591 for kpt in points[indices]:
592 for label, k in special_points.items():
593 if abs(kpt - k).sum() < eps:
594 break
595 else:
596 # No exact match. Try modulus 1:
597 for label, k in special_points.items():
598 if abs((kpt - k) % 1).sum() < eps:
599 break
600 else:
601 label = '?'
602 labels.append(label)
604 xcoords, ixcoords = indices_to_axis_coords(indices, points, cell)
605 return xcoords, ixcoords, labels
608def indices_to_axis_coords(indices, points, cell):
609 jump = False # marks a discontinuity in the path
610 xcoords = [0]
611 for i1, i2 in zip(indices[:-1], indices[1:]):
612 if not jump and i1 + 1 == i2:
613 length = 0
614 jump = True # we don't want two jumps in a row
615 else:
616 diff = points[i2] - points[i1]
617 length = np.linalg.norm(kpoint_convert(cell, skpts_kc=diff))
618 jump = False
619 xcoords.extend(np.linspace(0, length, i2 - i1 + 1)[1:] + xcoords[-1])
621 xcoords = np.array(xcoords)
622 return xcoords, xcoords[indices]
625special_paths = {
626 'cubic': 'GXMGRX,MR',
627 'fcc': 'GXWKGLUWLK,UX',
628 'bcc': 'GHNGPH,PN',
629 'tetragonal': 'GXMGZRAZXR,MA',
630 'orthorhombic': 'GXSYGZURTZ,YT,UX,SR',
631 'hexagonal': 'GMKGALHA,LM,KH',
632 'monoclinic': 'GYHCEM1AXH1,MDZ,YD',
633 'rhombohedral type 1': 'GLB1,BZGX,QFP1Z,LP',
634 'rhombohedral type 2': 'GPZQGFP1Q1LZ'}
637def get_special_points(cell, lattice=None, eps=2e-4):
638 """Return dict of special points.
640 The definitions are from a paper by Wahyu Setyawana and Stefano
641 Curtarolo::
643 https://doi.org/10.1016/j.commatsci.2010.05.010
645 cell: 3x3 ndarray
646 Unit cell.
647 lattice: str
648 Optionally check that the cell is one of the following: cubic, fcc,
649 bcc, orthorhombic, tetragonal, hexagonal or monoclinic.
650 eps: float
651 Tolerance for cell-check.
652 """
654 if isinstance(cell, str):
655 warnings.warn('Please call this function with cell as the first '
656 'argument')
657 lattice, cell = cell, lattice
659 cell = Cell.ascell(cell)
660 # We create the bandpath because we want to transform the kpoints too,
661 # from the canonical cell to the given one.
662 #
663 # Note that this function is missing a tolerance, epsilon.
664 path = cell.bandpath(npoints=0)
665 return path.special_points
668def monkhorst_pack_interpolate(path, values, icell, bz2ibz,
669 size, offset=(0, 0, 0), pad_width=2):
670 """Interpolate values from Monkhorst-Pack sampling.
672 `monkhorst_pack_interpolate` takes a set of `values`, for example
673 eigenvalues, that are resolved on a Monkhorst Pack k-point grid given by
674 `size` and `offset` and interpolates the values onto the k-points
675 in `path`.
677 Note
678 ----
679 For the interpolation to work, path has to lie inside the domain
680 that is spanned by the MP kpoint grid given by size and offset.
682 To try to ensure this we expand the domain slightly by adding additional
683 entries along the edges and sides of the domain with values determined by
684 wrapping the values to the opposite side of the domain. In this way we
685 assume that the function to be interpolated is a periodic function in
686 k-space. The padding width is determined by the `pad_width` parameter.
688 Parameters
689 ----------
690 path: (nk, 3) array-like
691 Desired path in units of reciprocal lattice vectors.
692 values: (nibz, ...) array-like
693 Values on Monkhorst-Pack grid.
694 icell: (3, 3) array-like
695 Reciprocal lattice vectors.
696 bz2ibz: (nbz,) array-like of int
697 Map from nbz points in BZ to nibz reduced points in IBZ.
698 size: (3,) array-like of int
699 Size of Monkhorst-Pack grid.
700 offset: (3,) array-like
701 Offset of Monkhorst-Pack grid.
702 pad_width: int
703 Padding width to aid interpolation
705 Returns
706 -------
707 (nbz,) array-like
708 *values* interpolated to *path*.
709 """
710 from scipy.interpolate import LinearNDInterpolator
712 path = (np.asarray(path) + 0.5) % 1 - 0.5
713 path = np.dot(path, icell)
715 # Fold out values from IBZ to BZ:
716 v = np.asarray(values)[bz2ibz]
717 v = v.reshape(tuple(size) + v.shape[1:])
719 # Create padded Monkhorst-Pack grid:
720 size = np.asarray(size)
721 i = (np.indices(size + 2 * pad_width)
722 .transpose((1, 2, 3, 0)).reshape((-1, 3)))
723 k = (i - pad_width + 0.5) / size - 0.5 + offset
724 k = np.dot(k, icell)
726 # Fill in boundary values:
727 V = np.pad(v, [(pad_width, pad_width)] * 3 +
728 [(0, 0)] * (v.ndim - 3), mode="wrap")
730 interpolate = LinearNDInterpolator(k, V.reshape((-1,) + V.shape[3:]))
731 interpolated_points = interpolate(path)
733 # NaN values indicate points outside interpolation domain, if fail
734 # try increasing padding
735 assert not np.isnan(interpolated_points).any(), \
736 "Points outside interpolation domain. Try increasing pad_width."
738 return interpolated_points
741# ChadiCohen k point grids. The k point grids are given in units of the
742# reciprocal unit cell. The variables are named after the following
743# convention: cc+'<Nkpoints>'+_+'shape'. For example an 18 k point
744# sq(3)xsq(3) is named 'cc18_sq3xsq3'.
746cc6_1x1 = np.array([
747 1, 1, 0, 1, 0, 0, 0, -1, 0, -1, -1, 0, -1, 0, 0,
748 0, 1, 0]).reshape((6, 3)) / 3.0
750cc12_2x3 = np.array([
751 3, 4, 0, 3, 10, 0, 6, 8, 0, 3, -2, 0, 6, -4, 0,
752 6, 2, 0, -3, 8, 0, -3, 2, 0, -3, -4, 0, -6, 4, 0, -6, -2, 0, -6,
753 -8, 0]).reshape((12, 3)) / 18.0
755cc18_sq3xsq3 = np.array([
756 2, 2, 0, 4, 4, 0, 8, 2, 0, 4, -2, 0, 8, -4,
757 0, 10, -2, 0, 10, -8, 0, 8, -10, 0, 2, -10, 0, 4, -8, 0, -2, -8,
758 0, 2, -4, 0, -4, -4, 0, -2, -2, 0, -4, 2, 0, -2, 4, 0, -8, 4, 0,
759 -4, 8, 0]).reshape((18, 3)) / 18.0
761cc18_1x1 = np.array([
762 2, 4, 0, 2, 10, 0, 4, 8, 0, 8, 4, 0, 8, 10, 0,
763 10, 8, 0, 2, -2, 0, 4, -4, 0, 4, 2, 0, -2, 8, 0, -2, 2, 0, -2, -4,
764 0, -4, 4, 0, -4, -2, 0, -4, -8, 0, -8, 2, 0, -8, -4, 0, -10, -2,
765 0]).reshape((18, 3)) / 18.0
767cc54_sq3xsq3 = np.array([
768 4, -10, 0, 6, -10, 0, 0, -8, 0, 2, -8, 0, 6,
769 -8, 0, 8, -8, 0, -4, -6, 0, -2, -6, 0, 2, -6, 0, 4, -6, 0, 8, -6,
770 0, 10, -6, 0, -6, -4, 0, -2, -4, 0, 0, -4, 0, 4, -4, 0, 6, -4, 0,
771 10, -4, 0, -6, -2, 0, -4, -2, 0, 0, -2, 0, 2, -2, 0, 6, -2, 0, 8,
772 -2, 0, -8, 0, 0, -4, 0, 0, -2, 0, 0, 2, 0, 0, 4, 0, 0, 8, 0, 0,
773 -8, 2, 0, -6, 2, 0, -2, 2, 0, 0, 2, 0, 4, 2, 0, 6, 2, 0, -10, 4,
774 0, -6, 4, 0, -4, 4, 0, 0, 4, 0, 2, 4, 0, 6, 4, 0, -10, 6, 0, -8,
775 6, 0, -4, 6, 0, -2, 6, 0, 2, 6, 0, 4, 6, 0, -8, 8, 0, -6, 8, 0,
776 -2, 8, 0, 0, 8, 0, -6, 10, 0, -4, 10, 0]).reshape((54, 3)) / 18.0
778cc54_1x1 = np.array([
779 2, 2, 0, 4, 4, 0, 8, 8, 0, 6, 8, 0, 4, 6, 0, 6,
780 10, 0, 4, 10, 0, 2, 6, 0, 2, 8, 0, 0, 2, 0, 0, 4, 0, 0, 8, 0, -2,
781 6, 0, -2, 4, 0, -4, 6, 0, -6, 4, 0, -4, 2, 0, -6, 2, 0, -2, 0, 0,
782 -4, 0, 0, -8, 0, 0, -8, -2, 0, -6, -2, 0, -10, -4, 0, -10, -6, 0,
783 -6, -4, 0, -8, -6, 0, -2, -2, 0, -4, -4, 0, -8, -8, 0, 4, -2, 0,
784 6, -2, 0, 6, -4, 0, 2, 0, 0, 4, 0, 0, 6, 2, 0, 6, 4, 0, 8, 6, 0,
785 8, 0, 0, 8, 2, 0, 10, 4, 0, 10, 6, 0, 2, -4, 0, 2, -6, 0, 4, -6,
786 0, 0, -2, 0, 0, -4, 0, -2, -6, 0, -4, -6, 0, -6, -8, 0, 0, -8, 0,
787 -2, -8, 0, -4, -10, 0, -6, -10, 0]).reshape((54, 3)) / 18.0
789cc162_sq3xsq3 = np.array([
790 -8, 16, 0, -10, 14, 0, -7, 14, 0, -4, 14,
791 0, -11, 13, 0, -8, 13, 0, -5, 13, 0, -2, 13, 0, -13, 11, 0, -10,
792 11, 0, -7, 11, 0, -4, 11, 0, -1, 11, 0, 2, 11, 0, -14, 10, 0, -11,
793 10, 0, -8, 10, 0, -5, 10, 0, -2, 10, 0, 1, 10, 0, 4, 10, 0, -16,
794 8, 0, -13, 8, 0, -10, 8, 0, -7, 8, 0, -4, 8, 0, -1, 8, 0, 2, 8, 0,
795 5, 8, 0, 8, 8, 0, -14, 7, 0, -11, 7, 0, -8, 7, 0, -5, 7, 0, -2, 7,
796 0, 1, 7, 0, 4, 7, 0, 7, 7, 0, 10, 7, 0, -13, 5, 0, -10, 5, 0, -7,
797 5, 0, -4, 5, 0, -1, 5, 0, 2, 5, 0, 5, 5, 0, 8, 5, 0, 11, 5, 0,
798 -14, 4, 0, -11, 4, 0, -8, 4, 0, -5, 4, 0, -2, 4, 0, 1, 4, 0, 4, 4,
799 0, 7, 4, 0, 10, 4, 0, -13, 2, 0, -10, 2, 0, -7, 2, 0, -4, 2, 0,
800 -1, 2, 0, 2, 2, 0, 5, 2, 0, 8, 2, 0, 11, 2, 0, -11, 1, 0, -8, 1,
801 0, -5, 1, 0, -2, 1, 0, 1, 1, 0, 4, 1, 0, 7, 1, 0, 10, 1, 0, 13, 1,
802 0, -10, -1, 0, -7, -1, 0, -4, -1, 0, -1, -1, 0, 2, -1, 0, 5, -1,
803 0, 8, -1, 0, 11, -1, 0, 14, -1, 0, -11, -2, 0, -8, -2, 0, -5, -2,
804 0, -2, -2, 0, 1, -2, 0, 4, -2, 0, 7, -2, 0, 10, -2, 0, 13, -2, 0,
805 -10, -4, 0, -7, -4, 0, -4, -4, 0, -1, -4, 0, 2, -4, 0, 5, -4, 0,
806 8, -4, 0, 11, -4, 0, 14, -4, 0, -8, -5, 0, -5, -5, 0, -2, -5, 0,
807 1, -5, 0, 4, -5, 0, 7, -5, 0, 10, -5, 0, 13, -5, 0, 16, -5, 0, -7,
808 -7, 0, -4, -7, 0, -1, -7, 0, 2, -7, 0, 5, -7, 0, 8, -7, 0, 11, -7,
809 0, 14, -7, 0, 17, -7, 0, -8, -8, 0, -5, -8, 0, -2, -8, 0, 1, -8,
810 0, 4, -8, 0, 7, -8, 0, 10, -8, 0, 13, -8, 0, 16, -8, 0, -7, -10,
811 0, -4, -10, 0, -1, -10, 0, 2, -10, 0, 5, -10, 0, 8, -10, 0, 11,
812 -10, 0, 14, -10, 0, 17, -10, 0, -5, -11, 0, -2, -11, 0, 1, -11, 0,
813 4, -11, 0, 7, -11, 0, 10, -11, 0, 13, -11, 0, 16, -11, 0, -1, -13,
814 0, 2, -13, 0, 5, -13, 0, 8, -13, 0, 11, -13, 0, 14, -13, 0, 1,
815 -14, 0, 4, -14, 0, 7, -14, 0, 10, -14, 0, 13, -14, 0, 5, -16, 0,
816 8, -16, 0, 11, -16, 0, 7, -17, 0, 10, -17, 0]).reshape((162, 3)) / 27.0
818cc162_1x1 = np.array([
819 -8, -16, 0, -10, -14, 0, -7, -14, 0, -4, -14,
820 0, -11, -13, 0, -8, -13, 0, -5, -13, 0, -2, -13, 0, -13, -11, 0,
821 -10, -11, 0, -7, -11, 0, -4, -11, 0, -1, -11, 0, 2, -11, 0, -14,
822 -10, 0, -11, -10, 0, -8, -10, 0, -5, -10, 0, -2, -10, 0, 1, -10,
823 0, 4, -10, 0, -16, -8, 0, -13, -8, 0, -10, -8, 0, -7, -8, 0, -4,
824 -8, 0, -1, -8, 0, 2, -8, 0, 5, -8, 0, 8, -8, 0, -14, -7, 0, -11,
825 -7, 0, -8, -7, 0, -5, -7, 0, -2, -7, 0, 1, -7, 0, 4, -7, 0, 7, -7,
826 0, 10, -7, 0, -13, -5, 0, -10, -5, 0, -7, -5, 0, -4, -5, 0, -1,
827 -5, 0, 2, -5, 0, 5, -5, 0, 8, -5, 0, 11, -5, 0, -14, -4, 0, -11,
828 -4, 0, -8, -4, 0, -5, -4, 0, -2, -4, 0, 1, -4, 0, 4, -4, 0, 7, -4,
829 0, 10, -4, 0, -13, -2, 0, -10, -2, 0, -7, -2, 0, -4, -2, 0, -1,
830 -2, 0, 2, -2, 0, 5, -2, 0, 8, -2, 0, 11, -2, 0, -11, -1, 0, -8,
831 -1, 0, -5, -1, 0, -2, -1, 0, 1, -1, 0, 4, -1, 0, 7, -1, 0, 10, -1,
832 0, 13, -1, 0, -10, 1, 0, -7, 1, 0, -4, 1, 0, -1, 1, 0, 2, 1, 0, 5,
833 1, 0, 8, 1, 0, 11, 1, 0, 14, 1, 0, -11, 2, 0, -8, 2, 0, -5, 2, 0,
834 -2, 2, 0, 1, 2, 0, 4, 2, 0, 7, 2, 0, 10, 2, 0, 13, 2, 0, -10, 4,
835 0, -7, 4, 0, -4, 4, 0, -1, 4, 0, 2, 4, 0, 5, 4, 0, 8, 4, 0, 11, 4,
836 0, 14, 4, 0, -8, 5, 0, -5, 5, 0, -2, 5, 0, 1, 5, 0, 4, 5, 0, 7, 5,
837 0, 10, 5, 0, 13, 5, 0, 16, 5, 0, -7, 7, 0, -4, 7, 0, -1, 7, 0, 2,
838 7, 0, 5, 7, 0, 8, 7, 0, 11, 7, 0, 14, 7, 0, 17, 7, 0, -8, 8, 0,
839 -5, 8, 0, -2, 8, 0, 1, 8, 0, 4, 8, 0, 7, 8, 0, 10, 8, 0, 13, 8, 0,
840 16, 8, 0, -7, 10, 0, -4, 10, 0, -1, 10, 0, 2, 10, 0, 5, 10, 0, 8,
841 10, 0, 11, 10, 0, 14, 10, 0, 17, 10, 0, -5, 11, 0, -2, 11, 0, 1,
842 11, 0, 4, 11, 0, 7, 11, 0, 10, 11, 0, 13, 11, 0, 16, 11, 0, -1,
843 13, 0, 2, 13, 0, 5, 13, 0, 8, 13, 0, 11, 13, 0, 14, 13, 0, 1, 14,
844 0, 4, 14, 0, 7, 14, 0, 10, 14, 0, 13, 14, 0, 5, 16, 0, 8, 16, 0,
845 11, 16, 0, 7, 17, 0, 10, 17, 0]).reshape((162, 3)) / 27.0
848# The following is a list of the critical points in the 1st Brillouin zone
849# for some typical crystal structures following the conventions of Setyawan
850# and Curtarolo [https://doi.org/10.1016/j.commatsci.2010.05.010].
851#
852# In units of the reciprocal basis vectors.
853#
854# See http://en.wikipedia.org/wiki/Brillouin_zone
855sc_special_points = {
856 'cubic': {'G': [0, 0, 0],
857 'M': [1 / 2, 1 / 2, 0],
858 'R': [1 / 2, 1 / 2, 1 / 2],
859 'X': [0, 1 / 2, 0]},
860 'fcc': {'G': [0, 0, 0],
861 'K': [3 / 8, 3 / 8, 3 / 4],
862 'L': [1 / 2, 1 / 2, 1 / 2],
863 'U': [5 / 8, 1 / 4, 5 / 8],
864 'W': [1 / 2, 1 / 4, 3 / 4],
865 'X': [1 / 2, 0, 1 / 2]},
866 'bcc': {'G': [0, 0, 0],
867 'H': [1 / 2, -1 / 2, 1 / 2],
868 'P': [1 / 4, 1 / 4, 1 / 4],
869 'N': [0, 0, 1 / 2]},
870 'tetragonal': {'G': [0, 0, 0],
871 'A': [1 / 2, 1 / 2, 1 / 2],
872 'M': [1 / 2, 1 / 2, 0],
873 'R': [0, 1 / 2, 1 / 2],
874 'X': [0, 1 / 2, 0],
875 'Z': [0, 0, 1 / 2]},
876 'orthorhombic': {'G': [0, 0, 0],
877 'R': [1 / 2, 1 / 2, 1 / 2],
878 'S': [1 / 2, 1 / 2, 0],
879 'T': [0, 1 / 2, 1 / 2],
880 'U': [1 / 2, 0, 1 / 2],
881 'X': [1 / 2, 0, 0],
882 'Y': [0, 1 / 2, 0],
883 'Z': [0, 0, 1 / 2]},
884 'hexagonal': {'G': [0, 0, 0],
885 'A': [0, 0, 1 / 2],
886 'H': [1 / 3, 1 / 3, 1 / 2],
887 'K': [1 / 3, 1 / 3, 0],
888 'L': [1 / 2, 0, 1 / 2],
889 'M': [1 / 2, 0, 0]}}
892# Old version of dictionary kept for backwards compatibility.
893# Not for ordinary use.
894ibz_points = {'cubic': {'Gamma': [0, 0, 0],
895 'X': [0, 0 / 2, 1 / 2],
896 'R': [1 / 2, 1 / 2, 1 / 2],
897 'M': [0 / 2, 1 / 2, 1 / 2]},
898 'fcc': {'Gamma': [0, 0, 0],
899 'X': [1 / 2, 0, 1 / 2],
900 'W': [1 / 2, 1 / 4, 3 / 4],
901 'K': [3 / 8, 3 / 8, 3 / 4],
902 'U': [5 / 8, 1 / 4, 5 / 8],
903 'L': [1 / 2, 1 / 2, 1 / 2]},
904 'bcc': {'Gamma': [0, 0, 0],
905 'H': [1 / 2, -1 / 2, 1 / 2],
906 'N': [0, 0, 1 / 2],
907 'P': [1 / 4, 1 / 4, 1 / 4]},
908 'hexagonal': {'Gamma': [0, 0, 0],
909 'M': [0, 1 / 2, 0],
910 'K': [-1 / 3, 1 / 3, 0],
911 'A': [0, 0, 1 / 2],
912 'L': [0, 1 / 2, 1 / 2],
913 'H': [-1 / 3, 1 / 3, 1 / 2]},
914 'tetragonal': {'Gamma': [0, 0, 0],
915 'X': [1 / 2, 0, 0],
916 'M': [1 / 2, 1 / 2, 0],
917 'Z': [0, 0, 1 / 2],
918 'R': [1 / 2, 0, 1 / 2],
919 'A': [1 / 2, 1 / 2, 1 / 2]},
920 'orthorhombic': {'Gamma': [0, 0, 0],
921 'R': [1 / 2, 1 / 2, 1 / 2],
922 'S': [1 / 2, 1 / 2, 0],
923 'T': [0, 1 / 2, 1 / 2],
924 'U': [1 / 2, 0, 1 / 2],
925 'X': [1 / 2, 0, 0],
926 'Y': [0, 1 / 2, 0],
927 'Z': [0, 0, 1 / 2]}}