Coverage for /builds/debichem-team/python-ase/ase/lattice/__init__.py : 96.07%
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
1from abc import abstractmethod, ABC
2import functools
3import warnings
4import numpy as np
5from typing import Dict, List
7from ase.cell import Cell
8from ase.build.bulk import bulk as newbulk
9from ase.dft.kpoints import parse_path_string, sc_special_points, BandPath
10from ase.utils import pbc2pbc
13@functools.wraps(newbulk)
14def bulk(*args, **kwargs):
15 warnings.warn('Use ase.build.bulk() instead', stacklevel=2)
16 return newbulk(*args, **kwargs)
19_degrees = np.pi / 180
22class BravaisLattice(ABC):
23 """Represent Bravais lattices and data related to the Brillouin zone.
25 There are 14 3D Bravais classes: CUB, FCC, BCC, ..., and TRI, and
26 five 2D classes.
28 Each class stores basic static information:
30 >>> from ase.lattice import FCC, MCL
31 >>> FCC.name
32 'FCC'
33 >>> FCC.longname
34 'face-centred cubic'
35 >>> FCC.pearson_symbol
36 'cF'
37 >>> MCL.parameters
38 ('a', 'b', 'c', 'alpha')
40 Each class can be instantiated with the specific lattice parameters
41 that apply to that lattice:
43 >>> MCL(3, 4, 5, 80)
44 MCL(a=3, b=4, c=5, alpha=80)
46 """
47 # These parameters can be set by the @bravais decorator for a subclass.
48 # (We could also use metaclasses to do this, but that's more abstract)
49 name = None # e.g. 'CUB', 'BCT', 'ORCF', ...
50 longname = None # e.g. 'cubic', 'body-centred tetragonal', ...
51 parameters = None # e.g. ('a', 'c')
52 variants = None # e.g. {'BCT1': <variant object>,
53 # 'BCT2': <variant object>}
54 ndim = None
56 def __init__(self, **kwargs):
57 p = {}
58 eps = kwargs.pop('eps', 2e-4)
59 for k, v in kwargs.items():
60 p[k] = float(v)
61 assert set(p) == set(self.parameters)
62 self._parameters = p
63 self._eps = eps
65 if len(self.variants) == 1:
66 # If there's only one it has the same name as the lattice type
67 self._variant = self.variants[self.name]
68 else:
69 name = self._variant_name(**self._parameters)
70 self._variant = self.variants[name]
72 @property
73 def variant(self) -> str:
74 """Return name of lattice variant.
76 >>> BCT(3, 5).variant
77 'BCT2'
78 """
79 return self._variant.name
81 def __getattr__(self, name: str):
82 if name in self._parameters:
83 return self._parameters[name]
84 return self.__getattribute__(name) # Raises error
86 def vars(self) -> Dict[str, float]:
87 """Get parameter names and values of this lattice as a dictionary."""
88 return dict(self._parameters)
90 def conventional(self) -> 'BravaisLattice':
91 """Get the conventional cell corresponding to this lattice."""
92 cls = bravais_lattices[self.conventional_cls]
93 return cls(**self._parameters)
95 def tocell(self) -> Cell:
96 """Return this lattice as a :class:`~ase.cell.Cell` object."""
97 cell = self._cell(**self._parameters)
98 return Cell(cell)
100 def get_transformation(self, cell, eps=1e-8):
101 # Get transformation matrix relating input cell to canonical cell
102 T = cell.dot(np.linalg.pinv(self.tocell()))
103 msg = 'This transformation changes the length/area/volume of the cell'
104 assert np.isclose(np.abs(np.linalg.det(T[:self.ndim,
105 :self.ndim])), 1,
106 atol=eps), msg
107 return T
109 def cellpar(self) -> np.ndarray:
110 """Get cell lengths and angles as array of length 6.
112 See :func:`ase.geometry.Cell.cellpar`."""
113 # (Just a brute-force implementation)
114 cell = self.tocell()
115 return cell.cellpar()
117 @property
118 def special_path(self) -> str:
119 """Get default special k-point path for this lattice as a string.
121 >>> BCT(3, 5).special_path
122 'GXYSGZS1NPY1Z,XP'
123 """
124 return self._variant.special_path
126 @property
127 def special_point_names(self) -> List[str]:
128 """Return all special point names as a list of strings.
130 >>> BCT(3, 5).special_point_names
131 ['G', 'N', 'P', 'S', 'S1', 'X', 'Y', 'Y1', 'Z']
132 """
133 labels = parse_path_string(self._variant.special_point_names)
134 assert len(labels) == 1 # list of lists
135 return labels[0]
137 def get_special_points_array(self) -> np.ndarray:
138 """Return all special points for this lattice as an array.
140 Ordering is consistent with special_point_names."""
141 if self._variant.special_points is not None:
142 # Fixed dictionary of special points
143 d = self.get_special_points()
144 labels = self.special_point_names
145 assert len(d) == len(labels)
146 points = np.empty((len(d), 3))
147 for i, label in enumerate(labels):
148 points[i] = d[label]
149 return points
151 # Special points depend on lattice parameters:
152 points = self._special_points(variant=self._variant,
153 **self._parameters)
154 assert len(points) == len(self.special_point_names)
155 return np.array(points)
157 def get_special_points(self) -> Dict[str, np.ndarray]:
158 """Return a dictionary of named special k-points for this lattice."""
159 if self._variant.special_points is not None:
160 return self._variant.special_points
162 labels = self.special_point_names
163 points = self.get_special_points_array()
165 return dict(zip(labels, points))
167 def plot_bz(self, path=None, special_points=None, **plotkwargs):
168 """Plot the reciprocal cell and default bandpath."""
169 # Create a generic bandpath (no interpolated kpoints):
170 bandpath = self.bandpath(path=path, special_points=special_points,
171 npoints=0)
172 return bandpath.plot(dimension=self.ndim, **plotkwargs)
174 def bandpath(self, path=None, npoints=None, special_points=None,
175 density=None, transformation=None) -> BandPath:
176 """Return a :class:`~ase.dft.kpoints.BandPath` for this lattice.
178 See :meth:`ase.cell.Cell.bandpath` for description of parameters.
180 >>> BCT(3, 5).bandpath()
181 BandPath(path='GXYSGZS1NPY1Z,XP', cell=[3x3], \
182special_points={GNPSS1XYY1Z}, kpts=[51x3])
184 .. note:: This produces the standard band path following AFlow
185 conventions. If your cell does not follow this convention,
186 you will need :meth:`ase.cell.Cell.bandpath` instead or
187 the kpoints may not correspond to your particular cell.
189 """
190 if special_points is None:
191 special_points = self.get_special_points()
193 if path is None:
194 path = self._variant.special_path
195 elif not isinstance(path, str):
196 from ase.dft.kpoints import resolve_custom_points
197 path, special_points = resolve_custom_points(path,
198 special_points,
199 self._eps)
201 cell = self.tocell()
202 if transformation is not None:
203 cell = transformation.dot(cell)
205 bandpath = BandPath(cell=cell, path=path,
206 special_points=special_points)
207 return bandpath.interpolate(npoints=npoints, density=density)
209 @abstractmethod
210 def _cell(self, **kwargs):
211 """Return a Cell object from this Bravais lattice.
213 Arguments are the dictionary of Bravais parameters."""
214 pass
216 def _special_points(self, **kwargs):
217 """Return the special point coordinates as an npoints x 3 sequence.
219 Subclasses typically return a nested list.
221 Ordering must be same as kpoint labels.
223 Arguments are the dictionary of Bravais parameters and the variant."""
224 raise NotImplementedError
226 def _variant_name(self, **kwargs):
227 """Return the name (e.g. 'ORCF3') of variant.
229 Arguments will be the dictionary of Bravais parameters."""
230 raise NotImplementedError
232 def __format__(self, spec):
233 tokens = []
234 if not spec:
235 spec = '.6g'
236 template = '{}={:%s}' % spec
238 for name in self.parameters:
239 value = self._parameters[name]
240 tokens.append(template.format(name, value))
241 return '{}({})'.format(self.name, ', '.join(tokens))
243 def __str__(self) -> str:
244 return self.__format__('')
246 def __repr__(self) -> str:
247 return self.__format__('.20g')
249 def description(self) -> str:
250 """Return complete description of lattice and Brillouin zone."""
251 points = self.get_special_points()
252 labels = self.special_point_names
254 coordstring = '\n'.join([' {:2s} {:7.4f} {:7.4f} {:7.4f}'
255 .format(label, *points[label])
256 for label in labels])
258 string = """\
259{repr}
260 {variant}
261 Special point coordinates:
262{special_points}
263""".format(repr=str(self),
264 variant=self._variant,
265 special_points=coordstring)
266 return string
268 @classmethod
269 def type_description(cls):
270 """Return complete description of this Bravais lattice type."""
271 desc = """\
272Lattice name: {name}
273 Long name: {longname}
274 Parameters: {parameters}
275""".format(**vars(cls))
277 chunks = [desc]
278 for name in cls.variant_names:
279 var = cls.variants[name]
280 txt = str(var)
281 lines = [' ' + L for L in txt.splitlines()]
282 lines.append('')
283 chunks.extend(lines)
285 return '\n'.join(chunks)
288class Variant:
289 variant_desc = """\
290Variant name: {name}
291 Special point names: {special_point_names}
292 Default path: {special_path}
293"""
295 def __init__(self, name, special_point_names, special_path,
296 special_points=None):
297 self.name = name
298 self.special_point_names = special_point_names
299 self.special_path = special_path
300 if special_points is not None:
301 special_points = dict(special_points)
302 for key, value in special_points.items():
303 special_points[key] = np.array(value)
304 self.special_points = special_points
305 # XXX Should make special_points available as a single array as well
306 # (easier to transform that way)
308 def __str__(self) -> str:
309 return self.variant_desc.format(**vars(self))
312bravais_names = []
313bravais_lattices = {}
314bravais_classes = {}
317def bravaisclass(longname, crystal_family, lattice_system, pearson_symbol,
318 parameters, variants, ndim=3):
319 """Decorator for Bravais lattice classes.
321 This sets a number of class variables and processes the information
322 about different variants into a list of Variant objects."""
324 def decorate(cls):
325 btype = cls.__name__
326 cls.name = btype
327 cls.longname = longname
328 cls.crystal_family = crystal_family
329 cls.lattice_system = lattice_system
330 cls.pearson_symbol = pearson_symbol
331 cls.parameters = tuple(parameters)
332 cls.variant_names = []
333 cls.variants = {}
334 cls.ndim = ndim
336 for [name, special_point_names, special_path,
337 special_points] in variants:
338 cls.variant_names.append(name)
339 cls.variants[name] = Variant(name, special_point_names,
340 special_path, special_points)
342 # Register in global list and dictionary
343 bravais_names.append(btype)
344 bravais_lattices[btype] = cls
345 bravais_classes[pearson_symbol] = cls
346 return cls
348 return decorate
351# Common mappings of primitive to conventional cells:
352_identity = np.identity(3, int)
353_fcc_map = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
354_bcc_map = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]])
357class UnconventionalLattice(ValueError):
358 pass
361class Cubic(BravaisLattice):
362 """Abstract class for cubic lattices."""
363 conventional_cls = 'CUB'
365 def __init__(self, a):
366 BravaisLattice.__init__(self, a=a)
369@bravaisclass('primitive cubic', 'cubic', 'cubic', 'cP', 'a',
370 [['CUB', 'GXRM', 'GXMGRX,MR', sc_special_points['cubic']]])
371class CUB(Cubic):
372 conventional_cellmap = _identity
374 def _cell(self, a):
375 return a * np.eye(3)
378@bravaisclass('face-centred cubic', 'cubic', 'cubic', 'cF', 'a',
379 [['FCC', 'GKLUWX', 'GXWKGLUWLK,UX', sc_special_points['fcc']]])
380class FCC(Cubic):
381 conventional_cellmap = _bcc_map
383 def _cell(self, a):
384 return 0.5 * np.array([[0., a, a], [a, 0, a], [a, a, 0]])
387@bravaisclass('body-centred cubic', 'cubic', 'cubic', 'cI', 'a',
388 [['BCC', 'GHPN', 'GHNGPH,PN', sc_special_points['bcc']]])
389class BCC(Cubic):
390 conventional_cellmap = _fcc_map
392 def _cell(self, a):
393 return 0.5 * np.array([[-a, a, a], [a, -a, a], [a, a, -a]])
396@bravaisclass('primitive tetragonal', 'tetragonal', 'tetragonal', 'tP', 'ac',
397 [['TET', 'GAMRXZ', 'GXMGZRAZ,XR,MA',
398 sc_special_points['tetragonal']]])
399class TET(BravaisLattice):
400 conventional_cls = 'TET'
401 conventional_cellmap = _identity
403 def __init__(self, a, c):
404 BravaisLattice.__init__(self, a=a, c=c)
406 def _cell(self, a, c):
407 return np.diag(np.array([a, a, c]))
410# XXX in BCT2 we use S for Sigma.
411# Also in other places I think
412@bravaisclass('body-centred tetragonal', 'tetragonal', 'tetragonal', 'tI',
413 'ac',
414 [['BCT1', 'GMNPXZZ1', 'GXMGZPNZ1M,XP', None],
415 ['BCT2', 'GNPSS1XYY1Z', 'GXYSGZS1NPY1Z,XP', None]])
416class BCT(BravaisLattice):
417 conventional_cls = 'TET'
418 conventional_cellmap = _fcc_map
420 def __init__(self, a, c):
421 BravaisLattice.__init__(self, a=a, c=c)
423 def _cell(self, a, c):
424 return 0.5 * np.array([[-a, a, c], [a, -a, c], [a, a, -c]])
426 def _variant_name(self, a, c):
427 return 'BCT1' if c < a else 'BCT2'
429 def _special_points(self, a, c, variant):
430 a2 = a * a
431 c2 = c * c
433 assert variant.name in self.variants
435 if variant.name == 'BCT1':
436 eta = .25 * (1 + c2 / a2)
437 points = [[0, 0, 0],
438 [-.5, .5, .5],
439 [0., .5, 0.],
440 [.25, .25, .25],
441 [0., 0., .5],
442 [eta, eta, -eta],
443 [-eta, 1 - eta, eta]]
444 else:
445 eta = .25 * (1 + a2 / c2) # Not same eta as BCT1!
446 zeta = 0.5 * a2 / c2
447 points = [[0., .0, 0.],
448 [0., .5, 0.],
449 [.25, .25, .25],
450 [-eta, eta, eta],
451 [eta, 1 - eta, -eta],
452 [0., 0., .5],
453 [-zeta, zeta, .5],
454 [.5, .5, -zeta],
455 [.5, .5, -.5]]
456 return points
459def check_orc(a, b, c):
460 if not a < b < c:
461 raise UnconventionalLattice('Expected a < b < c, got {}, {}, {}'
462 .format(a, b, c))
465class Orthorhombic(BravaisLattice):
466 """Abstract class for orthorhombic types."""
468 def __init__(self, a, b, c):
469 check_orc(a, b, c)
470 BravaisLattice.__init__(self, a=a, b=b, c=c)
473@bravaisclass('primitive orthorhombic', 'orthorhombic', 'orthorhombic', 'oP',
474 'abc',
475 [['ORC', 'GRSTUXYZ', 'GXSYGZURTZ,YT,UX,SR',
476 sc_special_points['orthorhombic']]])
477class ORC(Orthorhombic):
478 conventional_cls = 'ORC'
479 conventional_cellmap = _identity
481 def _cell(self, a, b, c):
482 return np.diag([a, b, c]).astype(float)
485@bravaisclass('face-centred orthorhombic', 'orthorhombic', 'orthorhombic',
486 'oF', 'abc',
487 [['ORCF1', 'GAA1LTXX1YZ', 'GYTZGXA1Y,TX1,XAZ,LG', None],
488 ['ORCF2', 'GCC1DD1LHH1XYZ', 'GYCDXGZD1HC,C1Z,XH1,HY,LG', None],
489 ['ORCF3', 'GAA1LTXX1YZ', 'GYTZGXA1Y,XAZ,LG', None]])
490class ORCF(Orthorhombic):
491 conventional_cls = 'ORC'
492 conventional_cellmap = _bcc_map
494 def _cell(self, a, b, c):
495 return 0.5 * np.array([[0, b, c], [a, 0, c], [a, b, 0]])
497 def _special_points(self, a, b, c, variant):
498 a2 = a * a
499 b2 = b * b
500 c2 = c * c
501 xminus = 0.25 * (1 + a2 / b2 - a2 / c2)
502 xplus = 0.25 * (1 + a2 / b2 + a2 / c2)
504 if variant.name == 'ORCF1' or variant.name == 'ORCF3':
505 zeta = xminus
506 eta = xplus
508 points = [[0, 0, 0],
509 [.5, .5 + zeta, zeta],
510 [.5, .5 - zeta, 1 - zeta],
511 [.5, .5, .5],
512 [1., .5, .5],
513 [0., eta, eta],
514 [1., 1 - eta, 1 - eta],
515 [.5, 0, .5],
516 [.5, .5, 0]]
517 else:
518 assert variant.name == 'ORCF2'
519 phi = 0.25 * (1 + c2 / b2 - c2 / a2)
520 delta = 0.25 * (1 + b2 / a2 - b2 / c2)
521 eta = xminus
523 points = [[0, 0, 0],
524 [.5, .5 - eta, 1 - eta],
525 [.5, .5 + eta, eta],
526 [.5 - delta, .5, 1 - delta],
527 [.5 + delta, .5, delta],
528 [.5, .5, .5],
529 [1 - phi, .5 - phi, .5],
530 [phi, .5 + phi, .5],
531 [0., .5, .5],
532 [.5, 0., .5],
533 [.5, .5, 0.]]
535 return points
537 def _variant_name(self, a, b, c):
538 diff = 1.0 / (a * a) - 1.0 / (b * b) - 1.0 / (c * c)
539 if abs(diff) < self._eps:
540 return 'ORCF3'
541 return 'ORCF1' if diff > 0 else 'ORCF2'
544@bravaisclass('body-centred orthorhombic', 'orthorhombic', 'orthorhombic',
545 'oI', 'abc',
546 [['ORCI', 'GLL1L2RSTWXX1YY1Z', 'GXLTWRX1ZGYSW,L1Y,Y1Z', None]])
547class ORCI(Orthorhombic):
548 conventional_cls = 'ORC'
549 conventional_cellmap = _fcc_map
551 def _cell(self, a, b, c):
552 return 0.5 * np.array([[-a, b, c], [a, -b, c], [a, b, -c]])
554 def _special_points(self, a, b, c, variant):
555 a2 = a**2
556 b2 = b**2
557 c2 = c**2
559 zeta = .25 * (1 + a2 / c2)
560 eta = .25 * (1 + b2 / c2)
561 delta = .25 * (b2 - a2) / c2
562 mu = .25 * (a2 + b2) / c2
564 points = [[0., 0., 0.],
565 [-mu, mu, .5 - delta],
566 [mu, -mu, .5 + delta],
567 [.5 - delta, .5 + delta, -mu],
568 [0, .5, 0],
569 [.5, 0, 0],
570 [0., 0., .5],
571 [.25, .25, .25],
572 [-zeta, zeta, zeta],
573 [zeta, 1 - zeta, -zeta],
574 [eta, -eta, eta],
575 [1 - eta, eta, -eta],
576 [.5, .5, -.5]]
577 return points
580@bravaisclass('base-centred orthorhombic', 'orthorhombic', 'orthorhombic',
581 'oC', 'abc',
582 [['ORCC', 'GAA1RSTXX1YZ', 'GXSRAZGYX1A1TY,ZT', None]])
583class ORCC(BravaisLattice):
584 conventional_cls = 'ORC'
585 conventional_cellmap = np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 1]])
587 def __init__(self, a, b, c):
588 # ORCC is the only ORCx lattice with a<b and not a<b<c
589 if a >= b:
590 raise UnconventionalLattice(f'Expected a < b, got a={a}, b={b}')
591 BravaisLattice.__init__(self, a=a, b=b, c=c)
593 def _cell(self, a, b, c):
594 return np.array([[0.5 * a, -0.5 * b, 0], [0.5 * a, 0.5 * b, 0],
595 [0, 0, c]])
597 def _special_points(self, a, b, c, variant):
598 zeta = .25 * (1 + a * a / (b * b))
599 points = [[0, 0, 0],
600 [zeta, zeta, .5],
601 [-zeta, 1 - zeta, .5],
602 [0, .5, .5],
603 [0, .5, 0],
604 [-.5, .5, .5],
605 [zeta, zeta, 0],
606 [-zeta, 1 - zeta, 0],
607 [-.5, .5, 0],
608 [0, 0, .5]]
609 return points
612@bravaisclass('primitive hexagonal', 'hexagonal', 'hexagonal', 'hP',
613 'ac',
614 [['HEX', 'GMKALH', 'GMKGALHA,LM,KH',
615 sc_special_points['hexagonal']]])
616class HEX(BravaisLattice):
617 conventional_cls = 'HEX'
618 conventional_cellmap = _identity
620 def __init__(self, a, c):
621 BravaisLattice.__init__(self, a=a, c=c)
623 def _cell(self, a, c):
624 x = 0.5 * np.sqrt(3)
625 return np.array([[0.5 * a, -x * a, 0], [0.5 * a, x * a, 0],
626 [0., 0., c]])
629@bravaisclass('primitive rhombohedral', 'hexagonal', 'rhombohedral', 'hR',
630 ('a', 'alpha'),
631 [['RHL1', 'GBB1FLL1PP1P2QXZ', 'GLB1,BZGX,QFP1Z,LP', None],
632 ['RHL2', 'GFLPP1QQ1Z', 'GPZQGFP1Q1LZ', None]])
633class RHL(BravaisLattice):
634 conventional_cls = 'RHL'
635 conventional_cellmap = _identity
637 def __init__(self, a, alpha):
638 if alpha >= 120:
639 raise UnconventionalLattice('Need alpha < 120 degrees, got {}'
640 .format(alpha))
641 BravaisLattice.__init__(self, a=a, alpha=alpha)
643 def _cell(self, a, alpha):
644 alpha *= np.pi / 180
645 acosa = a * np.cos(alpha)
646 acosa2 = a * np.cos(0.5 * alpha)
647 asina2 = a * np.sin(0.5 * alpha)
648 acosfrac = acosa / acosa2
649 xx = (1 - acosfrac**2)
650 assert xx > 0.0
651 return np.array([[acosa2, -asina2, 0], [acosa2, asina2, 0],
652 [a * acosfrac, 0, a * xx**0.5]])
654 def _variant_name(self, a, alpha):
655 return 'RHL1' if alpha < 90 else 'RHL2'
657 def _special_points(self, a, alpha, variant):
658 if variant.name == 'RHL1':
659 cosa = np.cos(alpha * _degrees)
660 eta = (1 + 4 * cosa) / (2 + 4 * cosa)
661 nu = .75 - 0.5 * eta
662 points = [[0, 0, 0],
663 [eta, .5, 1 - eta],
664 [.5, 1 - eta, eta - 1],
665 [.5, .5, 0],
666 [.5, 0, 0],
667 [0, 0, -.5],
668 [eta, nu, nu],
669 [1 - nu, 1 - nu, 1 - eta],
670 [nu, nu, eta - 1],
671 [1 - nu, nu, 0],
672 [nu, 0, -nu],
673 [.5, .5, .5]]
674 else:
675 eta = 1 / (2 * np.tan(alpha * _degrees / 2)**2)
676 nu = .75 - 0.5 * eta
677 points = [[0, 0, 0],
678 [.5, -.5, 0],
679 [.5, 0, 0],
680 [1 - nu, -nu, 1 - nu],
681 [nu, nu - 1, nu - 1],
682 [eta, eta, eta],
683 [1 - eta, -eta, -eta],
684 [.5, -.5, .5]]
685 return points
688def check_mcl(a, b, c, alpha):
689 if not (b <= c and alpha < 90):
690 raise UnconventionalLattice('Expected b <= c, alpha < 90; '
691 'got a={}, b={}, c={}, alpha={}'
692 .format(a, b, c, alpha))
695@bravaisclass('primitive monoclinic', 'monoclinic', 'monoclinic', 'mP',
696 ('a', 'b', 'c', 'alpha'),
697 [['MCL', 'GACDD1EHH1H2MM1M2XYY1Z', 'GYHCEM1AXH1,MDZ,YD', None]])
698class MCL(BravaisLattice):
699 conventional_cls = 'MCL'
700 conventional_cellmap = _identity
702 def __init__(self, a, b, c, alpha):
703 check_mcl(a, b, c, alpha)
704 BravaisLattice.__init__(self, a=a, b=b, c=c, alpha=alpha)
706 def _cell(self, a, b, c, alpha):
707 alpha *= _degrees
708 return np.array([[a, 0, 0], [0, b, 0],
709 [0, c * np.cos(alpha), c * np.sin(alpha)]])
711 def _special_points(self, a, b, c, alpha, variant):
712 cosa = np.cos(alpha * _degrees)
713 eta = (1 - b * cosa / c) / (2 * np.sin(alpha * _degrees)**2)
714 nu = .5 - eta * c * cosa / b
716 points = [[0, 0, 0],
717 [.5, .5, 0],
718 [0, .5, .5],
719 [.5, 0, .5],
720 [.5, 0, -.5],
721 [.5, .5, .5],
722 [0, eta, 1 - nu],
723 [0, 1 - eta, nu],
724 [0, eta, -nu],
725 [.5, eta, 1 - nu],
726 [.5, 1 - eta, nu],
727 [.5, eta, -nu],
728 [0, .5, 0],
729 [0, 0, .5],
730 [0, 0, -.5],
731 [.5, 0, 0]]
732 return points
734 def _variant_name(self, a, b, c, alpha):
735 check_mcl(a, b, c, alpha)
736 return 'MCL'
739@bravaisclass('base-centred monoclinic', 'monoclinic', 'monoclinic', 'mC',
740 ('a', 'b', 'c', 'alpha'),
741 [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z',
742 'GYFLI,I1ZF1,YX1,XGN,MG', None],
743 ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z',
744 'GYFLI,I1ZF1,NGM', None],
745 ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',
746 'GYFHZIF1,H1Y1XGN,MG', None],
747 ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',
748 'GYFHZI,H1Y1XGN,MG', None],
749 ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z',
750 'GYFLI,I1ZHF1,H1Y1XGN,MG', None]])
751class MCLC(BravaisLattice):
752 conventional_cls = 'MCL'
753 conventional_cellmap = np.array([[1, -1, 0], [1, 1, 0], [0, 0, 1]])
755 def __init__(self, a, b, c, alpha):
756 check_mcl(a, b, c, alpha)
757 BravaisLattice.__init__(self, a=a, b=b, c=c, alpha=alpha)
759 def _cell(self, a, b, c, alpha):
760 alpha *= np.pi / 180
761 return np.array([[0.5 * a, 0.5 * b, 0], [-0.5 * a, 0.5 * b, 0],
762 [0, c * np.cos(alpha), c * np.sin(alpha)]])
764 def _variant_name(self, a, b, c, alpha):
765 # from ase.geometry.cell import mclc
766 # okay, this is a bit hacky
768 # We need the same parameters here as when determining the points.
769 # Right now we just repeat the code:
770 check_mcl(a, b, c, alpha)
772 a2 = a * a
773 b2 = b * b
774 cosa = np.cos(alpha * _degrees)
775 sina = np.sin(alpha * _degrees)
776 sina2 = sina**2
778 cell = self.tocell()
779 lengths_angles = Cell(cell.reciprocal()).cellpar()
781 kgamma = lengths_angles[-1]
783 eps = self._eps
784 # We should not compare angles in degrees versus lengths with
785 # the same precision.
786 if abs(kgamma - 90) < eps:
787 variant = 2
788 elif kgamma > 90:
789 variant = 1
790 elif kgamma < 90:
791 num = b * cosa / c + b2 * sina2 / a2
792 if abs(num - 1) < eps:
793 variant = 4
794 elif num < 1:
795 variant = 3
796 else:
797 variant = 5
798 variant = 'MCLC' + str(variant)
799 return variant
801 def _special_points(self, a, b, c, alpha, variant):
802 variant = int(variant.name[-1])
804 a2 = a * a
805 b2 = b * b
806 # c2 = c * c
807 cosa = np.cos(alpha * _degrees)
808 sina = np.sin(alpha * _degrees)
809 sina2 = sina**2
811 if variant == 1 or variant == 2:
812 zeta = (2 - b * cosa / c) / (4 * sina2)
813 eta = 0.5 + 2 * zeta * c * cosa / b
814 psi = .75 - a2 / (4 * b2 * sina * sina)
815 phi = psi + (.75 - psi) * b * cosa / c
817 points = [[0, 0, 0],
818 [.5, 0, 0],
819 [0, -.5, 0],
820 [1 - zeta, 1 - zeta, 1 - eta],
821 [zeta, zeta, eta],
822 [-zeta, -zeta, 1 - eta],
823 [1 - zeta, -zeta, 1 - eta],
824 [phi, 1 - phi, .5],
825 [1 - phi, phi - 1, .5],
826 [.5, .5, .5],
827 [.5, 0, .5],
828 [1 - psi, psi - 1, 0],
829 [psi, 1 - psi, 0],
830 [psi - 1, -psi, 0],
831 [.5, .5, 0],
832 [-.5, -.5, 0],
833 [0, 0, .5]]
834 elif variant == 3 or variant == 4:
835 mu = .25 * (1 + b2 / a2)
836 delta = b * c * cosa / (2 * a2)
837 zeta = mu - 0.25 + (1 - b * cosa / c) / (4 * sina2)
838 eta = 0.5 + 2 * zeta * c * cosa / b
839 phi = 1 + zeta - 2 * mu
840 psi = eta - 2 * delta
842 points = [[0, 0, 0],
843 [1 - phi, 1 - phi, 1 - psi],
844 [phi, phi - 1, psi],
845 [1 - phi, -phi, 1 - psi],
846 [zeta, zeta, eta],
847 [1 - zeta, -zeta, 1 - eta],
848 [-zeta, -zeta, 1 - eta],
849 [.5, -.5, .5],
850 [.5, 0, .5],
851 [.5, 0, 0],
852 [0, -.5, 0],
853 [.5, -.5, 0],
854 [mu, mu, delta],
855 [1 - mu, -mu, -delta],
856 [-mu, -mu, -delta],
857 [mu, mu - 1, delta],
858 [0, 0, .5]]
859 elif variant == 5:
860 zeta = .25 * (b2 / a2 + (1 - b * cosa / c) / sina2)
861 eta = 0.5 + 2 * zeta * c * cosa / b
862 mu = .5 * eta + b2 / (4 * a2) - b * c * cosa / (2 * a2)
863 nu = 2 * mu - zeta
864 omega = (4 * nu - 1 - b2 * sina2 / a2) * c / (2 * b * cosa)
865 delta = zeta * c * cosa / b + omega / 2 - .25
866 rho = 1 - zeta * a2 / b2
868 points = [[0, 0, 0],
869 [nu, nu, omega],
870 [1 - nu, 1 - nu, 1 - omega],
871 [nu, nu - 1, omega],
872 [zeta, zeta, eta],
873 [1 - zeta, -zeta, 1 - eta],
874 [-zeta, -zeta, 1 - eta],
875 [rho, 1 - rho, .5],
876 [1 - rho, rho - 1, .5],
877 [.5, .5, .5],
878 [.5, 0, .5],
879 [.5, 0, 0],
880 [0, -.5, 0],
881 [.5, -.5, 0],
882 [mu, mu, delta],
883 [1 - mu, -mu, -delta],
884 [-mu, -mu, -delta],
885 [mu, mu - 1, delta],
886 [0, 0, .5]]
888 return points
891tri_angles_explanation = """\
892Angles kalpha, kbeta and kgamma of TRI lattice must be 1) all greater \
893than 90 degrees with kgamma being the smallest, or 2) all smaller than \
89490 with kgamma being the largest, or 3) kgamma=90 being the \
895smallest of the three, or 4) kgamma=90 being the largest of the three. \
896Angles of reciprocal lattice are kalpha={}, kbeta={}, kgamma={}. \
897If you don't care, please use Cell.fromcellpar() instead."""
899# XXX labels, paths, are all the same.
902@bravaisclass('primitive triclinic', 'triclinic', 'triclinic', 'aP',
903 ('a', 'b', 'c', 'alpha', 'beta', 'gamma'),
904 [['TRI1a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
905 ['TRI2a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
906 ['TRI1b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
907 ['TRI2b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None]])
908class TRI(BravaisLattice):
909 conventional_cls = 'TRI'
910 conventional_cellmap = _identity
912 def __init__(self, a, b, c, alpha, beta, gamma):
913 BravaisLattice.__init__(self, a=a, b=b, c=c, alpha=alpha, beta=beta,
914 gamma=gamma)
916 def _cell(self, a, b, c, alpha, beta, gamma):
917 alpha, beta, gamma = np.array([alpha, beta, gamma])
918 singamma = np.sin(gamma * _degrees)
919 cosgamma = np.cos(gamma * _degrees)
920 cosbeta = np.cos(beta * _degrees)
921 cosalpha = np.cos(alpha * _degrees)
922 a3x = c * cosbeta
923 a3y = c / singamma * (cosalpha - cosbeta * cosgamma)
924 a3z = c / singamma * np.sqrt(singamma**2 - cosalpha**2 - cosbeta**2
925 + 2 * cosalpha * cosbeta * cosgamma)
926 return np.array([[a, 0, 0], [b * cosgamma, b * singamma, 0],
927 [a3x, a3y, a3z]])
929 def _variant_name(self, a, b, c, alpha, beta, gamma):
930 cell = Cell.new([a, b, c, alpha, beta, gamma])
931 icellpar = Cell(cell.reciprocal()).cellpar()
932 kangles = kalpha, kbeta, kgamma = icellpar[3:]
934 def raise_unconventional():
935 raise UnconventionalLattice(tri_angles_explanation
936 .format(*kangles))
938 eps = self._eps
939 if abs(kgamma - 90) < eps:
940 if kalpha > 90 and kbeta > 90:
941 var = '2a'
942 elif kalpha < 90 and kbeta < 90:
943 var = '2b'
944 else:
945 # Is this possible? Maybe due to epsilon
946 raise_unconventional()
947 elif all(kangles > 90):
948 if kgamma > min(kangles):
949 raise_unconventional()
950 var = '1a'
951 elif all(kangles < 90): # and kgamma > max(kalpha, kbeta):
952 if kgamma < max(kangles):
953 raise_unconventional()
954 var = '1b'
955 else:
956 raise_unconventional()
958 return 'TRI' + var
960 def _special_points(self, a, b, c, alpha, beta, gamma, variant):
961 # (None of the points actually depend on any parameters)
962 # (We should store the points openly on the variant objects)
963 if variant.name == 'TRI1a' or variant.name == 'TRI2a':
964 points = [[0., 0., 0.],
965 [.5, .5, 0],
966 [0, .5, .5],
967 [.5, 0, .5],
968 [.5, .5, .5],
969 [.5, 0, 0],
970 [0, .5, 0],
971 [0, 0, .5]]
972 else:
973 points = [[0, 0, 0],
974 [.5, -.5, 0],
975 [0, 0, .5],
976 [-.5, -.5, .5],
977 [0, -.5, .5],
978 [0, -0.5, 0],
979 [.5, 0, 0],
980 [-.5, 0, .5]]
982 return points
985def get_subset_points(names, points):
986 newpoints = {}
987 for name in names:
988 newpoints[name] = points[name]
990 return newpoints
993@bravaisclass('primitive oblique', 'monoclinic', None, 'mp',
994 ('a', 'b', 'alpha'), [['OBL', 'GYHCH1X', 'GYHCH1XG', None]],
995 ndim=2)
996class OBL(BravaisLattice):
997 def __init__(self, a, b, alpha, **kwargs):
998 BravaisLattice.__init__(self, a=a, b=b, alpha=alpha, **kwargs)
1000 def _cell(self, a, b, alpha):
1001 cosa = np.cos(alpha * _degrees)
1002 sina = np.sin(alpha * _degrees)
1004 return np.array([[a, 0, 0],
1005 [b * cosa, b * sina, 0],
1006 [0., 0., 0.]])
1008 def _special_points(self, a, b, alpha, variant):
1009 # XXX Check me
1010 if alpha > 90:
1011 _alpha = 180 - alpha
1012 a, b = b, a
1013 else:
1014 _alpha = alpha
1016 cosa = np.cos(_alpha * _degrees)
1017 eta = (1 - a * cosa / b) / (2 * np.sin(_alpha * _degrees)**2)
1018 nu = .5 - eta * b * cosa / a
1020 points = [[0, 0, 0],
1021 [0, 0.5, 0],
1022 [eta, 1 - nu, 0],
1023 [.5, .5, 0],
1024 [1 - eta, nu, 0],
1025 [.5, 0, 0]]
1027 if alpha > 90:
1028 # Then we map the special points back
1029 op = np.array([[0, 1, 0],
1030 [-1, 0, 0],
1031 [0, 0, 1]])
1032 points = np.dot(points, op.T)
1034 return points
1037@bravaisclass('primitive hexagonal', 'hexagonal', None, 'hp', 'a',
1038 [['HEX2D', 'GMK', 'GMKG',
1039 get_subset_points('GMK',
1040 sc_special_points['hexagonal'])]],
1041 ndim=2)
1042class HEX2D(BravaisLattice):
1043 def __init__(self, a, **kwargs):
1044 BravaisLattice.__init__(self, a=a, **kwargs)
1046 def _cell(self, a):
1047 x = 0.5 * np.sqrt(3)
1048 return np.array([[a, 0, 0],
1049 [-0.5 * a, x * a, 0],
1050 [0., 0., 0.]])
1053@bravaisclass('primitive rectangular', 'orthorhombic', None, 'op', 'ab',
1054 [['RECT', 'GXSY', 'GXSYGS',
1055 get_subset_points('GXSY',
1056 sc_special_points['orthorhombic'])]],
1057 ndim=2)
1058class RECT(BravaisLattice):
1059 def __init__(self, a, b, **kwargs):
1060 BravaisLattice.__init__(self, a=a, b=b, **kwargs)
1062 def _cell(self, a, b):
1063 return np.array([[a, 0, 0],
1064 [0, b, 0],
1065 [0, 0, 0.]])
1068@bravaisclass('centred rectangular', 'orthorhombic', None, 'oc',
1069 ('a', 'alpha'), [['CRECT', 'GXA1Y', 'GXA1YG', None]], ndim=2)
1070class CRECT(BravaisLattice):
1071 def __init__(self, a, alpha, **kwargs):
1072 BravaisLattice.__init__(self, a=a, alpha=alpha, **kwargs)
1074 def _cell(self, a, alpha):
1075 x = np.cos(alpha * _degrees)
1076 y = np.sin(alpha * _degrees)
1077 return np.array([[a, 0, 0],
1078 [a * x, a * y, 0],
1079 [0, 0, 0.]])
1081 def _special_points(self, a, alpha, variant):
1082 if alpha > 90:
1083 _alpha = 180 - alpha
1084 else:
1085 _alpha = alpha
1086 sina2 = np.sin(_alpha / 2 * _degrees)**2
1087 sina = np.sin(_alpha * _degrees)**2
1088 eta = sina2 / sina
1089 cosa = np.cos(_alpha * _degrees)
1090 xi = eta * cosa
1092 points = [[0, 0, 0],
1093 [eta, - eta, 0],
1094 [0.5 + xi, 0.5 - xi, 0],
1095 [0.5, 0.5, 0]]
1097 if alpha > 90:
1098 # Then we map the special points back
1099 op = np.array([[0, 1, 0],
1100 [-1, 0, 0],
1101 [0, 0, 1]])
1102 points = np.dot(points, op.T)
1103 return points
1106@bravaisclass('primitive square', 'tetragonal', None, 'tp', ('a',),
1107 [['SQR', 'GMX', 'MGXM',
1108 get_subset_points('GMX', sc_special_points['tetragonal'])]],
1109 ndim=2)
1110class SQR(BravaisLattice):
1111 def __init__(self, a, **kwargs):
1112 BravaisLattice.__init__(self, a=a, **kwargs)
1114 def _cell(self, a):
1115 return np.array([[a, 0, 0],
1116 [0, a, 0],
1117 [0, 0, 0.]])
1120@bravaisclass('primitive line', 'line', None, '?', ('a',),
1121 [['LINE', 'GX', 'GX', {'G': [0, 0, 0], 'X': [0.5, 0, 0]}]],
1122 ndim=1)
1123class LINE(BravaisLattice):
1124 def __init__(self, a, **kwargs):
1125 BravaisLattice.__init__(self, a=a, **kwargs)
1127 def _cell(self, a):
1128 return np.array([[a, 0.0, 0.0],
1129 [0.0, 0.0, 0.0],
1130 [0.0, 0.0, 0.0]])
1133def celldiff(cell1, cell2):
1134 """Return a unitless measure of the difference between two cells."""
1135 cell1 = Cell.ascell(cell1).complete()
1136 cell2 = Cell.ascell(cell2).complete()
1137 v1v2 = cell1.volume * cell2.volume
1138 if v1v2 == 0:
1139 raise ZeroDivisionError('Cell volumes are zero')
1140 scale = v1v2**(-1. / 3.) # --> 1/Ang^2
1141 x1 = cell1 @ cell1.T
1142 x2 = cell2 @ cell2.T
1143 dev = scale * np.abs(x2 - x1).max()
1144 return dev
1147def get_lattice_from_canonical_cell(cell, eps=2e-4):
1148 """Return a Bravais lattice representing the given cell.
1150 This works only for cells that are derived from the standard form
1151 (as generated by lat.tocell()) or rotations thereof.
1153 If the given cell does not resemble the known form of a Bravais
1154 lattice, raise RuntimeError."""
1155 return LatticeChecker(cell, eps).match()
1158def identify_lattice(cell, eps=2e-4, *, pbc=True):
1159 """Find Bravais lattice representing this cell.
1161 Returns Bravais lattice object representing the cell along with
1162 an operation that, applied to the cell, yields the same lengths
1163 and angles as the Bravais lattice object."""
1165 pbc = cell.any(1) & pbc2pbc(pbc)
1166 npbc = sum(pbc)
1168 if npbc == 1:
1169 i = np.argmax(pbc) # index of periodic axis
1170 a = cell[i, i]
1171 if a < 0 or cell[i, [i - 1, i - 2]].any():
1172 raise ValueError('Not a 1-d cell ASE can handle: {cell}.'
1173 .format(cell=cell))
1174 if i == 0:
1175 op = np.eye(3)
1176 elif i == 1:
1177 op = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
1178 else:
1179 op = np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]])
1180 return LINE(a), op
1182 if npbc == 2:
1183 lat, op = get_2d_bravais_lattice(cell, eps, pbc=pbc)
1184 return lat, op
1186 if npbc != 3:
1187 raise ValueError('System must be periodic either '
1188 'along all three axes, '
1189 'along two first axes or, '
1190 'along the third axis. '
1191 'Got pbc={}'.format(pbc))
1193 from ase.geometry.bravais_type_engine import niggli_op_table
1195 if cell.rank < 3:
1196 raise ValueError('Expected 3 linearly independent cell vectors')
1197 rcell, reduction_op = cell.niggli_reduce(eps=eps)
1199 # We tabulate the cell's Niggli-mapped versions so we don't need to
1200 # redo any work when the same Niggli-operation appears multiple times
1201 # in the table:
1202 memory = {}
1204 # We loop through the most symmetric kinds (CUB etc.) and return
1205 # the first one we find:
1206 for latname in LatticeChecker.check_order:
1207 # There may be multiple Niggli operations that produce valid
1208 # lattices, at least for MCL. In that case we will pick the
1209 # one whose angle is closest to 90, but it means we cannot
1210 # just return the first one we find so we must remember then:
1211 matching_lattices = []
1213 for op_key in niggli_op_table[latname]:
1214 checker_and_op = memory.get(op_key)
1215 if checker_and_op is None:
1216 normalization_op = np.array(op_key).reshape(3, 3)
1217 candidate = Cell(np.linalg.inv(normalization_op.T) @ rcell)
1218 checker = LatticeChecker(candidate, eps=eps)
1219 memory[op_key] = (checker, normalization_op)
1220 else:
1221 checker, normalization_op = checker_and_op
1223 lat = checker.query(latname)
1224 if lat is not None:
1225 op = normalization_op @ np.linalg.inv(reduction_op)
1226 matching_lattices.append((lat, op))
1228 # Among any matching lattices, return the one with lowest
1229 # orthogonality defect:
1230 best = None
1231 best_defect = np.inf
1232 for lat, op in matching_lattices:
1233 cell = lat.tocell()
1234 lengths = cell.lengths()
1235 defect = np.prod(lengths) / cell.volume
1236 if defect < best_defect:
1237 best = lat, op
1238 best_defect = defect
1240 if best is not None:
1241 return best
1244class LatticeChecker:
1245 # The check order is slightly different than elsewhere listed order
1246 # as we need to check HEX/RHL before the ORCx family.
1247 check_order = ['CUB', 'FCC', 'BCC', 'TET', 'BCT', 'HEX', 'RHL',
1248 'ORC', 'ORCF', 'ORCI', 'ORCC', 'MCL', 'MCLC', 'TRI']
1250 def __init__(self, cell, eps=2e-4):
1251 """Generate Bravais lattices that look (or not) like the given cell.
1253 The cell must be reduced to canonical form, i.e., it must
1254 be possible to produce a cell with the same lengths and angles
1255 by directly through one of the Bravais lattice classes.
1257 Generally for internal use (this module).
1259 For each of the 14 Bravais lattices, this object can produce
1260 a lattice object which represents the same cell, or None if
1261 the tolerance eps is not met."""
1262 self.cell = cell
1263 self.eps = eps
1265 self.cellpar = cell.cellpar()
1266 self.lengths = self.A, self.B, self.C = self.cellpar[:3]
1267 self.angles = self.cellpar[3:]
1269 # Use a 'neutral' length for checking cubic lattices
1270 self.A0 = self.lengths.mean()
1272 # Vector of the diagonal and then off-diagonal dot products:
1273 # [a1 · a1, a2 · a2, a3 · a3, a2 · a3, a3 · a1, a1 · a2]
1274 self.prods = (cell @ cell.T).flat[[0, 4, 8, 5, 2, 1]]
1276 def _check(self, latcls, *args):
1277 if any(arg <= 0 for arg in args):
1278 return None
1279 try:
1280 lat = latcls(*args)
1281 except UnconventionalLattice:
1282 return None
1284 newcell = lat.tocell()
1285 err = celldiff(self.cell, newcell)
1286 if err < self.eps:
1287 return lat
1289 def match(self):
1290 """Match cell against all lattices, returning most symmetric match.
1292 Returns the lattice object. Raises RuntimeError on failure."""
1293 for name in self.check_order:
1294 lat = self.query(name)
1295 if lat:
1296 return lat
1297 else:
1298 raise RuntimeError('Could not find lattice type for cell '
1299 'with lengths and angles {}'
1300 .format(self.cell.cellpar().tolist()))
1302 def query(self, latname):
1303 """Match cell against named Bravais lattice.
1305 Return lattice object on success, None on failure."""
1306 meth = getattr(self, latname)
1307 lat = meth()
1308 return lat
1310 def CUB(self):
1311 # These methods (CUB, FCC, ...) all return a lattice object if
1312 # it matches, else None.
1313 return self._check(CUB, self.A0)
1315 def FCC(self):
1316 return self._check(FCC, np.sqrt(2) * self.A0)
1318 def BCC(self):
1319 return self._check(BCC, 2.0 * self.A0 / np.sqrt(3))
1321 def TET(self):
1322 return self._check(TET, self.A, self.C)
1324 def _bct_orci_lengths(self):
1325 # Coordinate-system independent relation for BCT and ORCI
1326 # standard cells:
1327 # a1 · a1 + a2 · a3 == a² / 2
1328 # a2 · a2 + a3 · a1 == a² / 2 (BCT)
1329 # == b² / 2 (ORCI)
1330 # a3 · a3 + a1 · a2 == c² / 2
1331 # We use these to get a, b, and c in those cases.
1332 prods = self.prods
1333 lengthsqr = 2.0 * (prods[:3] + prods[3:])
1334 if any(lengthsqr < 0):
1335 return None
1336 return np.sqrt(lengthsqr)
1338 def BCT(self):
1339 lengths = self._bct_orci_lengths()
1340 if lengths is None:
1341 return None
1342 return self._check(BCT, lengths[0], lengths[2])
1344 def HEX(self):
1345 return self._check(HEX, self.A, self.C)
1347 def RHL(self):
1348 return self._check(RHL, self.A, self.angles[0])
1350 def ORC(self):
1351 return self._check(ORC, *self.lengths)
1353 def ORCF(self):
1354 # ORCF standard cell:
1355 # a2 · a3 = a²/4
1356 # a3 · a1 = b²/4
1357 # a1 · a2 = c²/4
1358 prods = self.prods
1359 if all(prods[3:] > 0):
1360 orcf_abc = 2 * np.sqrt(prods[3:])
1361 return self._check(ORCF, *orcf_abc)
1363 def ORCI(self):
1364 lengths = self._bct_orci_lengths()
1365 if lengths is None:
1366 return None
1367 return self._check(ORCI, *lengths)
1369 def _orcc_ab(self):
1370 # ORCC: a1 · a1 + a2 · a3 = a²/2
1371 # a2 · a2 - a2 · a3 = b²/2
1372 prods = self.prods
1373 orcc_sqr_ab = np.empty(2)
1374 orcc_sqr_ab[0] = 2.0 * (prods[0] + prods[5])
1375 orcc_sqr_ab[1] = 2.0 * (prods[1] - prods[5])
1376 if all(orcc_sqr_ab > 0):
1377 return np.sqrt(orcc_sqr_ab)
1379 def ORCC(self):
1380 orcc_lengths_ab = self._orcc_ab()
1381 if orcc_lengths_ab is None:
1382 return None
1383 return self._check(ORCC, *orcc_lengths_ab, self.C)
1385 def MCL(self):
1386 return self._check(MCL, *self.lengths, self.angles[0])
1388 def MCLC(self):
1389 # MCLC is similar to ORCC:
1390 orcc_ab = self._orcc_ab()
1391 if orcc_ab is None:
1392 return None
1394 prods = self.prods
1395 C = self.C
1396 mclc_a, mclc_b = orcc_ab[::-1] # a, b reversed wrt. ORCC
1397 mclc_cosa = 2.0 * prods[3] / (mclc_b * C)
1398 if -1 < mclc_cosa < 1:
1399 mclc_alpha = np.arccos(mclc_cosa) * 180 / np.pi
1400 if mclc_b > C:
1401 # XXX Temporary fix for certain otherwise
1402 # unrecognizable lattices.
1403 #
1404 # This error could happen if the input lattice maps to
1405 # something just outside the domain of conventional
1406 # lattices (less than the tolerance). Our solution is to
1407 # propose a nearby conventional lattice instead, which
1408 # will then be accepted if it's close enough.
1409 mclc_b = 0.5 * (mclc_b + C)
1410 C = mclc_b
1411 return self._check(MCLC, mclc_a, mclc_b, C, mclc_alpha)
1413 def TRI(self):
1414 return self._check(TRI, *self.cellpar)
1417class UnsupportedLattice(ValueError):
1418 pass
1421def get_2d_bravais_lattice(origcell, eps=2e-4, *, pbc=True):
1423 pbc = origcell.any(1) & pbc2pbc(pbc)
1424 if list(pbc) != [1, 1, 0]:
1425 raise UnsupportedLattice('Can only get 2D Bravais lattice of cell with '
1426 'pbc==[1, 1, 0]; but we have {}'.format(pbc))
1428 nonperiodic = pbc.argmin()
1429 # Start with op = I
1430 ops = [np.eye(3)]
1431 for i in range(-1, 1):
1432 for j in range(-1, 1):
1433 op = [[1, j],
1434 [i, 1]]
1435 if np.abs(np.linalg.det(op)) > 1e-5:
1436 # Only touch periodic dirs:
1437 op = np.insert(op, nonperiodic, [0, 0], 0)
1438 op = np.insert(op, nonperiodic, ~pbc, 1)
1439 ops.append(np.array(op))
1441 def allclose(a, b):
1442 return np.allclose(a, b, atol=eps)
1444 symrank = 0
1445 for op in ops:
1446 cell = Cell(op.dot(origcell))
1447 cellpar = cell.cellpar()
1448 angles = cellpar[3:]
1449 # Find a, b and gamma
1450 gamma = angles[~pbc][0]
1451 a, b = cellpar[:3][pbc]
1453 anglesm90 = np.abs(angles - 90)
1454 # Maximum one angle different from 90 deg in 2d please
1455 if np.sum(anglesm90 > eps) > 1:
1456 continue
1458 all_lengths_equal = abs(a - b) < eps
1460 if all_lengths_equal:
1461 if allclose(gamma, 90):
1462 lat = SQR(a)
1463 rank = 5
1464 elif allclose(gamma, 120):
1465 lat = HEX2D(a)
1466 rank = 4
1467 else:
1468 lat = CRECT(a, gamma)
1469 rank = 3
1470 else:
1471 if allclose(gamma, 90):
1472 lat = RECT(a, b)
1473 rank = 2
1474 else:
1475 lat = OBL(a, b, gamma)
1476 rank = 1
1478 op = lat.get_transformation(origcell, eps=eps)
1479 if not allclose(np.dot(op, lat.tocell())[pbc][:, pbc],
1480 origcell.array[pbc][:, pbc]):
1481 msg = ('Cannot recognize cell at all somehow! {}, {}, {}'.
1482 format(a, b, gamma))
1483 raise RuntimeError(msg)
1484 if rank > symrank:
1485 finalop = op
1486 symrank = rank
1487 finallat = lat
1489 return finallat, finalop.T
1492def all_variants(include_blunt_angles=True):
1493 """For testing and examples; yield all variants of all lattices."""
1494 a, b, c = 3., 4., 5.
1495 alpha = 55.0
1496 yield CUB(a)
1497 yield FCC(a)
1498 yield BCC(a)
1499 yield TET(a, c)
1500 bct1 = BCT(2 * a, c)
1501 bct2 = BCT(a, c)
1502 assert bct1.variant == 'BCT1'
1503 assert bct2.variant == 'BCT2'
1505 yield bct1
1506 yield bct2
1508 yield ORC(a, b, c)
1510 a0 = np.sqrt(1.0 / (1 / b**2 + 1 / c**2))
1511 orcf1 = ORCF(0.5 * a0, b, c)
1512 orcf2 = ORCF(1.2 * a0, b, c)
1513 orcf3 = ORCF(a0, b, c)
1514 assert orcf1.variant == 'ORCF1'
1515 assert orcf2.variant == 'ORCF2'
1516 assert orcf3.variant == 'ORCF3'
1517 yield orcf1
1518 yield orcf2
1519 yield orcf3
1521 yield ORCI(a, b, c)
1522 yield ORCC(a, b, c)
1524 yield HEX(a, c)
1526 rhl1 = RHL(a, alpha=55.0)
1527 assert rhl1.variant == 'RHL1'
1528 yield rhl1
1530 rhl2 = RHL(a, alpha=105.0)
1531 assert rhl2.variant == 'RHL2'
1532 yield rhl2
1534 # With these lengths, alpha < 65 (or so) would result in a lattice that
1535 # could also be represented with alpha > 65, which is more conventional.
1536 yield MCL(a, b, c, alpha=70.0)
1538 mclc1 = MCLC(a, b, c, 80)
1539 assert mclc1.variant == 'MCLC1'
1540 yield mclc1
1541 # mclc2 has same special points as mclc1
1543 mclc3 = MCLC(1.8 * a, b, c * 2, 80)
1544 assert mclc3.variant == 'MCLC3'
1545 yield mclc3
1546 # mclc4 has same special points as mclc3
1548 # XXX We should add MCLC2 and MCLC4 as well.
1550 mclc5 = MCLC(b, b, 1.1 * b, 70)
1551 assert mclc5.variant == 'MCLC5'
1552 yield mclc5
1554 def get_tri(kcellpar):
1555 # We build the TRI lattices from cellpars of reciprocal cell
1556 icell = Cell.fromcellpar(kcellpar)
1557 cellpar = Cell(4 * icell.reciprocal()).cellpar()
1558 return TRI(*cellpar)
1560 tri1a = get_tri([1., 1.2, 1.4, 120., 110., 100.])
1561 assert tri1a.variant == 'TRI1a'
1562 yield tri1a
1564 tri1b = get_tri([1., 1.2, 1.4, 50., 60., 70.])
1565 assert tri1b.variant == 'TRI1b'
1566 yield tri1b
1568 tri2a = get_tri([1., 1.2, 1.4, 120., 110., 90.])
1569 assert tri2a.variant == 'TRI2a'
1570 yield tri2a
1571 tri2b = get_tri([1., 1.2, 1.4, 50., 60., 90.])
1572 assert tri2b.variant == 'TRI2b'
1573 yield tri2b
1575 yield OBL(a, b, alpha=alpha)
1576 yield RECT(a, b)
1577 yield CRECT(a, alpha=alpha)
1578 yield HEX2D(a)
1579 yield SQR(a)
1580 yield LINE(a)
1582 if include_blunt_angles:
1583 beta = 110
1584 yield OBL(a, b, alpha=beta)
1585 yield CRECT(a, alpha=beta)