Coverage for /builds/ase/ase/ase/lattice/__init__.py : 96.77%
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 cellpar(self) -> np.ndarray:
101 """Get cell lengths and angles as array of length 6.
103 See :func:`ase.geometry.Cell.cellpar`."""
104 # (Just a brute-force implementation)
105 cell = self.tocell()
106 return cell.cellpar()
108 @property
109 def special_path(self) -> str:
110 """Get default special k-point path for this lattice as a string.
112 >>> BCT(3, 5).special_path
113 'GXYSGZS1NPY1Z,XP'
114 """
115 return self._variant.special_path
117 @property
118 def special_point_names(self) -> List[str]:
119 """Return all special point names as a list of strings.
121 >>> BCT(3, 5).special_point_names
122 ['G', 'N', 'P', 'S', 'S1', 'X', 'Y', 'Y1', 'Z']
123 """
124 labels = parse_path_string(self._variant.special_point_names)
125 assert len(labels) == 1 # list of lists
126 return labels[0]
128 def get_special_points_array(self) -> np.ndarray:
129 """Return all special points for this lattice as an array.
131 Ordering is consistent with special_point_names."""
132 if self._variant.special_points is not None:
133 # Fixed dictionary of special points
134 d = self.get_special_points()
135 labels = self.special_point_names
136 assert len(d) == len(labels)
137 points = np.empty((len(d), 3))
138 for i, label in enumerate(labels):
139 points[i] = d[label]
140 return points
142 # Special points depend on lattice parameters:
143 points = self._special_points(variant=self._variant,
144 **self._parameters)
145 assert len(points) == len(self.special_point_names)
146 return np.array(points)
148 def get_special_points(self) -> Dict[str, np.ndarray]:
149 """Return a dictionary of named special k-points for this lattice."""
150 if self._variant.special_points is not None:
151 return self._variant.special_points
153 labels = self.special_point_names
154 points = self.get_special_points_array()
156 return dict(zip(labels, points))
158 def plot_bz(self, path=None, special_points=None, **plotkwargs):
159 """Plot the reciprocal cell and default bandpath."""
160 # Create a generic bandpath (no interpolated kpoints):
161 bandpath = self.bandpath(path=path, special_points=special_points,
162 npoints=0)
163 return bandpath.plot(dimension=self.ndim, **plotkwargs)
165 def bandpath(self, path=None, npoints=None, special_points=None,
166 density=None) -> BandPath:
167 """Return a :class:`~ase.dft.kpoints.BandPath` for this lattice.
169 See :meth:`ase.cell.Cell.bandpath` for description of parameters.
171 >>> BCT(3, 5).bandpath()
172 BandPath(path='GXYSGZS1NPY1Z,XP', cell=[3x3], \
173special_points={GNPSS1XYY1Z}, kpts=[51x3])
175 .. note:: This produces the standard band path following AFlow
176 conventions. If your cell does not follow this convention,
177 you will need :meth:`ase.cell.Cell.bandpath` instead or
178 the kpoints may not correspond to your particular cell.
180 """
181 if special_points is None:
182 special_points = self.get_special_points()
184 if path is None:
185 path = self._variant.special_path
186 elif not isinstance(path, str):
187 from ase.dft.kpoints import resolve_custom_points
188 path, special_points = resolve_custom_points(path,
189 special_points,
190 self._eps)
192 cell = self.tocell()
193 bandpath = BandPath(cell=cell, path=path,
194 special_points=special_points)
195 return bandpath.interpolate(npoints=npoints, density=density)
197 @abstractmethod
198 def _cell(self, **kwargs):
199 """Return a Cell object from this Bravais lattice.
201 Arguments are the dictionary of Bravais parameters."""
203 def _special_points(self, **kwargs):
204 """Return the special point coordinates as an npoints x 3 sequence.
206 Subclasses typically return a nested list.
208 Ordering must be same as kpoint labels.
210 Arguments are the dictionary of Bravais parameters and the variant."""
211 raise NotImplementedError
213 def _variant_name(self, **kwargs):
214 """Return the name (e.g. 'ORCF3') of variant.
216 Arguments will be the dictionary of Bravais parameters."""
217 raise NotImplementedError
219 def __format__(self, spec):
220 tokens = []
221 if not spec:
222 spec = '.6g'
223 template = '{}={:%s}' % spec
225 for name in self.parameters:
226 value = self._parameters[name]
227 tokens.append(template.format(name, value))
228 return '{}({})'.format(self.name, ', '.join(tokens))
230 def __str__(self) -> str:
231 return self.__format__('')
233 def __repr__(self) -> str:
234 return self.__format__('.20g')
236 def description(self) -> str:
237 """Return complete description of lattice and Brillouin zone."""
238 points = self.get_special_points()
239 labels = self.special_point_names
241 coordstring = '\n'.join([' {:2s} {:7.4f} {:7.4f} {:7.4f}'
242 .format(label, *points[label])
243 for label in labels])
245 string = """\
246{repr}
247 {variant}
248 Special point coordinates:
249{special_points}
250""".format(repr=str(self),
251 variant=self._variant,
252 special_points=coordstring)
253 return string
255 @classmethod
256 def type_description(cls):
257 """Return complete description of this Bravais lattice type."""
258 desc = """\
259Lattice name: {name}
260 Long name: {longname}
261 Parameters: {parameters}
262""".format(**vars(cls))
264 chunks = [desc]
265 for name in cls.variant_names:
266 var = cls.variants[name]
267 txt = str(var)
268 lines = [' ' + L for L in txt.splitlines()]
269 lines.append('')
270 chunks.extend(lines)
272 return '\n'.join(chunks)
275class Variant:
276 variant_desc = """\
277Variant name: {name}
278 Special point names: {special_point_names}
279 Default path: {special_path}
280"""
282 def __init__(self, name, special_point_names, special_path,
283 special_points=None):
284 self.name = name
285 self.special_point_names = special_point_names
286 self.special_path = special_path
287 if special_points is not None:
288 special_points = dict(special_points)
289 for key, value in special_points.items():
290 special_points[key] = np.array(value)
291 self.special_points = special_points
292 # XXX Should make special_points available as a single array as well
293 # (easier to transform that way)
295 def __str__(self) -> str:
296 return self.variant_desc.format(**vars(self))
299bravais_names = []
300bravais_lattices = {}
301bravais_classes = {}
304def bravaisclass(longname, crystal_family, lattice_system, pearson_symbol,
305 parameters, variants, ndim=3):
306 """Decorator for Bravais lattice classes.
308 This sets a number of class variables and processes the information
309 about different variants into a list of Variant objects."""
311 def decorate(cls):
312 btype = cls.__name__
313 cls.name = btype
314 cls.longname = longname
315 cls.crystal_family = crystal_family
316 cls.lattice_system = lattice_system
317 cls.pearson_symbol = pearson_symbol
318 cls.parameters = tuple(parameters)
319 cls.variant_names = []
320 cls.variants = {}
321 cls.ndim = ndim
323 for [name, special_point_names, special_path,
324 special_points] in variants:
325 cls.variant_names.append(name)
326 cls.variants[name] = Variant(name, special_point_names,
327 special_path, special_points)
329 # Register in global list and dictionary
330 bravais_names.append(btype)
331 bravais_lattices[btype] = cls
332 bravais_classes[pearson_symbol] = cls
333 return cls
335 return decorate
338# Common mappings of primitive to conventional cells:
339_identity = np.identity(3, int)
340_fcc_map = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
341_bcc_map = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]])
344class UnconventionalLattice(ValueError):
345 pass
348class Cubic(BravaisLattice):
349 """Abstract class for cubic lattices."""
350 conventional_cls = 'CUB'
352 def __init__(self, a):
353 super().__init__(a=a)
356@bravaisclass('primitive cubic', 'cubic', 'cubic', 'cP', 'a',
357 [['CUB', 'GXRM', 'GXMGRX,MR', sc_special_points['cubic']]])
358class CUB(Cubic):
359 conventional_cellmap = _identity
361 def _cell(self, a):
362 return a * np.eye(3)
365@bravaisclass('face-centred cubic', 'cubic', 'cubic', 'cF', 'a',
366 [['FCC', 'GKLUWX', 'GXWKGLUWLK,UX', sc_special_points['fcc']]])
367class FCC(Cubic):
368 conventional_cellmap = _bcc_map
370 def _cell(self, a):
371 return 0.5 * np.array([[0., a, a], [a, 0, a], [a, a, 0]])
374@bravaisclass('body-centred cubic', 'cubic', 'cubic', 'cI', 'a',
375 [['BCC', 'GHPN', 'GHNGPH,PN', sc_special_points['bcc']]])
376class BCC(Cubic):
377 conventional_cellmap = _fcc_map
379 def _cell(self, a):
380 return 0.5 * np.array([[-a, a, a], [a, -a, a], [a, a, -a]])
383@bravaisclass('primitive tetragonal', 'tetragonal', 'tetragonal', 'tP', 'ac',
384 [['TET', 'GAMRXZ', 'GXMGZRAZ,XR,MA',
385 sc_special_points['tetragonal']]])
386class TET(BravaisLattice):
387 conventional_cls = 'TET'
388 conventional_cellmap = _identity
390 def __init__(self, a, c):
391 super().__init__(a=a, c=c)
393 def _cell(self, a, c):
394 return np.diag(np.array([a, a, c]))
397# XXX in BCT2 we use S for Sigma.
398# Also in other places I think
399@bravaisclass('body-centred tetragonal', 'tetragonal', 'tetragonal', 'tI',
400 'ac',
401 [['BCT1', 'GMNPXZZ1', 'GXMGZPNZ1M,XP', None],
402 ['BCT2', 'GNPSS1XYY1Z', 'GXYSGZS1NPY1Z,XP', None]])
403class BCT(BravaisLattice):
404 conventional_cls = 'TET'
405 conventional_cellmap = _fcc_map
407 def __init__(self, a, c):
408 super().__init__(a=a, c=c)
410 def _cell(self, a, c):
411 return 0.5 * np.array([[-a, a, c], [a, -a, c], [a, a, -c]])
413 def _variant_name(self, a, c):
414 return 'BCT1' if c < a else 'BCT2'
416 def _special_points(self, a, c, variant):
417 a2 = a * a
418 c2 = c * c
420 assert variant.name in self.variants
422 if variant.name == 'BCT1':
423 eta = .25 * (1 + c2 / a2)
424 points = [[0, 0, 0],
425 [-.5, .5, .5],
426 [0., .5, 0.],
427 [.25, .25, .25],
428 [0., 0., .5],
429 [eta, eta, -eta],
430 [-eta, 1 - eta, eta]]
431 else:
432 eta = .25 * (1 + a2 / c2) # Not same eta as BCT1!
433 zeta = 0.5 * a2 / c2
434 points = [[0., .0, 0.],
435 [0., .5, 0.],
436 [.25, .25, .25],
437 [-eta, eta, eta],
438 [eta, 1 - eta, -eta],
439 [0., 0., .5],
440 [-zeta, zeta, .5],
441 [.5, .5, -zeta],
442 [.5, .5, -.5]]
443 return points
446def check_orc(a, b, c):
447 if not a < b < c:
448 raise UnconventionalLattice('Expected a < b < c, got {}, {}, {}'
449 .format(a, b, c))
452class Orthorhombic(BravaisLattice):
453 """Abstract class for orthorhombic types."""
455 def __init__(self, a, b, c):
456 check_orc(a, b, c)
457 super().__init__(a=a, b=b, c=c)
460@bravaisclass('primitive orthorhombic', 'orthorhombic', 'orthorhombic', 'oP',
461 'abc',
462 [['ORC', 'GRSTUXYZ', 'GXSYGZURTZ,YT,UX,SR',
463 sc_special_points['orthorhombic']]])
464class ORC(Orthorhombic):
465 conventional_cls = 'ORC'
466 conventional_cellmap = _identity
468 def _cell(self, a, b, c):
469 return np.diag([a, b, c]).astype(float)
472@bravaisclass('face-centred orthorhombic', 'orthorhombic', 'orthorhombic',
473 'oF', 'abc',
474 [['ORCF1', 'GAA1LTXX1YZ', 'GYTZGXA1Y,TX1,XAZ,LG', None],
475 ['ORCF2', 'GCC1DD1LHH1XYZ', 'GYCDXGZD1HC,C1Z,XH1,HY,LG', None],
476 ['ORCF3', 'GAA1LTXX1YZ', 'GYTZGXA1Y,XAZ,LG', None]])
477class ORCF(Orthorhombic):
478 conventional_cls = 'ORC'
479 conventional_cellmap = _bcc_map
481 def _cell(self, a, b, c):
482 return 0.5 * np.array([[0, b, c], [a, 0, c], [a, b, 0]])
484 def _special_points(self, a, b, c, variant):
485 a2 = a * a
486 b2 = b * b
487 c2 = c * c
488 xminus = 0.25 * (1 + a2 / b2 - a2 / c2)
489 xplus = 0.25 * (1 + a2 / b2 + a2 / c2)
491 if variant.name == 'ORCF1' or variant.name == 'ORCF3':
492 zeta = xminus
493 eta = xplus
495 points = [[0, 0, 0],
496 [.5, .5 + zeta, zeta],
497 [.5, .5 - zeta, 1 - zeta],
498 [.5, .5, .5],
499 [1., .5, .5],
500 [0., eta, eta],
501 [1., 1 - eta, 1 - eta],
502 [.5, 0, .5],
503 [.5, .5, 0]]
504 else:
505 assert variant.name == 'ORCF2'
506 phi = 0.25 * (1 + c2 / b2 - c2 / a2)
507 delta = 0.25 * (1 + b2 / a2 - b2 / c2)
508 eta = xminus
510 points = [[0, 0, 0],
511 [.5, .5 - eta, 1 - eta],
512 [.5, .5 + eta, eta],
513 [.5 - delta, .5, 1 - delta],
514 [.5 + delta, .5, delta],
515 [.5, .5, .5],
516 [1 - phi, .5 - phi, .5],
517 [phi, .5 + phi, .5],
518 [0., .5, .5],
519 [.5, 0., .5],
520 [.5, .5, 0.]]
522 return points
524 def _variant_name(self, a, b, c):
525 diff = 1.0 / (a * a) - 1.0 / (b * b) - 1.0 / (c * c)
526 if abs(diff) < self._eps:
527 return 'ORCF3'
528 return 'ORCF1' if diff > 0 else 'ORCF2'
531@bravaisclass('body-centred orthorhombic', 'orthorhombic', 'orthorhombic',
532 'oI', 'abc',
533 [['ORCI', 'GLL1L2RSTWXX1YY1Z', 'GXLTWRX1ZGYSW,L1Y,Y1Z', None]])
534class ORCI(Orthorhombic):
535 conventional_cls = 'ORC'
536 conventional_cellmap = _fcc_map
538 def _cell(self, a, b, c):
539 return 0.5 * np.array([[-a, b, c], [a, -b, c], [a, b, -c]])
541 def _special_points(self, a, b, c, variant):
542 a2 = a**2
543 b2 = b**2
544 c2 = c**2
546 zeta = .25 * (1 + a2 / c2)
547 eta = .25 * (1 + b2 / c2)
548 delta = .25 * (b2 - a2) / c2
549 mu = .25 * (a2 + b2) / c2
551 points = [[0., 0., 0.],
552 [-mu, mu, .5 - delta],
553 [mu, -mu, .5 + delta],
554 [.5 - delta, .5 + delta, -mu],
555 [0, .5, 0],
556 [.5, 0, 0],
557 [0., 0., .5],
558 [.25, .25, .25],
559 [-zeta, zeta, zeta],
560 [zeta, 1 - zeta, -zeta],
561 [eta, -eta, eta],
562 [1 - eta, eta, -eta],
563 [.5, .5, -.5]]
564 return points
567@bravaisclass('base-centred orthorhombic', 'orthorhombic', 'orthorhombic',
568 'oC', 'abc',
569 [['ORCC', 'GAA1RSTXX1YZ', 'GXSRAZGYX1A1TY,ZT', None]])
570class ORCC(BravaisLattice):
571 conventional_cls = 'ORC'
572 conventional_cellmap = np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 1]])
574 def __init__(self, a, b, c):
575 # ORCC is the only ORCx lattice with a<b and not a<b<c
576 if a >= b:
577 raise UnconventionalLattice(f'Expected a < b, got a={a}, b={b}')
578 super().__init__(a=a, b=b, c=c)
580 def _cell(self, a, b, c):
581 return np.array([[0.5 * a, -0.5 * b, 0], [0.5 * a, 0.5 * b, 0],
582 [0, 0, c]])
584 def _special_points(self, a, b, c, variant):
585 zeta = .25 * (1 + a * a / (b * b))
586 points = [[0, 0, 0],
587 [zeta, zeta, .5],
588 [-zeta, 1 - zeta, .5],
589 [0, .5, .5],
590 [0, .5, 0],
591 [-.5, .5, .5],
592 [zeta, zeta, 0],
593 [-zeta, 1 - zeta, 0],
594 [-.5, .5, 0],
595 [0, 0, .5]]
596 return points
599@bravaisclass('primitive hexagonal', 'hexagonal', 'hexagonal', 'hP',
600 'ac',
601 [['HEX', 'GMKALH', 'GMKGALHA,LM,KH',
602 sc_special_points['hexagonal']]])
603class HEX(BravaisLattice):
604 conventional_cls = 'HEX'
605 conventional_cellmap = _identity
607 def __init__(self, a, c):
608 super().__init__(a=a, c=c)
610 def _cell(self, a, c):
611 x = 0.5 * np.sqrt(3)
612 return np.array([[0.5 * a, -x * a, 0], [0.5 * a, x * a, 0],
613 [0., 0., c]])
616@bravaisclass('primitive rhombohedral', 'hexagonal', 'rhombohedral', 'hR',
617 ('a', 'alpha'),
618 [['RHL1', 'GBB1FLL1PP1P2QXZ', 'GLB1,BZGX,QFP1Z,LP', None],
619 ['RHL2', 'GFLPP1QQ1Z', 'GPZQGFP1Q1LZ', None]])
620class RHL(BravaisLattice):
621 conventional_cls = 'RHL'
622 conventional_cellmap = _identity
624 def __init__(self, a, alpha):
625 if alpha >= 120:
626 raise UnconventionalLattice('Need alpha < 120 degrees, got {}'
627 .format(alpha))
628 super().__init__(a=a, alpha=alpha)
630 def _cell(self, a, alpha):
631 alpha *= np.pi / 180
632 acosa = a * np.cos(alpha)
633 acosa2 = a * np.cos(0.5 * alpha)
634 asina2 = a * np.sin(0.5 * alpha)
635 acosfrac = acosa / acosa2
636 xx = (1 - acosfrac**2)
637 assert xx > 0.0
638 return np.array([[acosa2, -asina2, 0], [acosa2, asina2, 0],
639 [a * acosfrac, 0, a * xx**0.5]])
641 def _variant_name(self, a, alpha):
642 return 'RHL1' if alpha < 90 else 'RHL2'
644 def _special_points(self, a, alpha, variant):
645 if variant.name == 'RHL1':
646 cosa = np.cos(alpha * _degrees)
647 eta = (1 + 4 * cosa) / (2 + 4 * cosa)
648 nu = .75 - 0.5 * eta
649 points = [[0, 0, 0],
650 [eta, .5, 1 - eta],
651 [.5, 1 - eta, eta - 1],
652 [.5, .5, 0],
653 [.5, 0, 0],
654 [0, 0, -.5],
655 [eta, nu, nu],
656 [1 - nu, 1 - nu, 1 - eta],
657 [nu, nu, eta - 1],
658 [1 - nu, nu, 0],
659 [nu, 0, -nu],
660 [.5, .5, .5]]
661 else:
662 eta = 1 / (2 * np.tan(alpha * _degrees / 2)**2)
663 nu = .75 - 0.5 * eta
664 points = [[0, 0, 0],
665 [.5, -.5, 0],
666 [.5, 0, 0],
667 [1 - nu, -nu, 1 - nu],
668 [nu, nu - 1, nu - 1],
669 [eta, eta, eta],
670 [1 - eta, -eta, -eta],
671 [.5, -.5, .5]]
672 return points
675def check_mcl(a, b, c, alpha):
676 if not (b <= c and alpha < 90):
677 raise UnconventionalLattice('Expected b <= c, alpha < 90; '
678 'got a={}, b={}, c={}, alpha={}'
679 .format(a, b, c, alpha))
682@bravaisclass('primitive monoclinic', 'monoclinic', 'monoclinic', 'mP',
683 ('a', 'b', 'c', 'alpha'),
684 [['MCL', 'GACDD1EHH1H2MM1M2XYY1Z', 'GYHCEM1AXH1,MDZ,YD', None]])
685class MCL(BravaisLattice):
686 conventional_cls = 'MCL'
687 conventional_cellmap = _identity
689 def __init__(self, a, b, c, alpha):
690 check_mcl(a, b, c, alpha)
691 super().__init__(a=a, b=b, c=c, alpha=alpha)
693 def _cell(self, a, b, c, alpha):
694 alpha *= _degrees
695 return np.array([[a, 0, 0], [0, b, 0],
696 [0, c * np.cos(alpha), c * np.sin(alpha)]])
698 def _special_points(self, a, b, c, alpha, variant):
699 cosa = np.cos(alpha * _degrees)
700 eta = (1 - b * cosa / c) / (2 * np.sin(alpha * _degrees)**2)
701 nu = .5 - eta * c * cosa / b
703 points = [[0, 0, 0],
704 [.5, .5, 0],
705 [0, .5, .5],
706 [.5, 0, .5],
707 [.5, 0, -.5],
708 [.5, .5, .5],
709 [0, eta, 1 - nu],
710 [0, 1 - eta, nu],
711 [0, eta, -nu],
712 [.5, eta, 1 - nu],
713 [.5, 1 - eta, nu],
714 [.5, eta, -nu],
715 [0, .5, 0],
716 [0, 0, .5],
717 [0, 0, -.5],
718 [.5, 0, 0]]
719 return points
721 def _variant_name(self, a, b, c, alpha):
722 check_mcl(a, b, c, alpha)
723 return 'MCL'
726@bravaisclass('base-centred monoclinic', 'monoclinic', 'monoclinic', 'mC',
727 ('a', 'b', 'c', 'alpha'),
728 [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z',
729 'GYFLI,I1ZF1,YX1,XGN,MG', None],
730 ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z',
731 'GYFLI,I1ZF1,NGM', None],
732 ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',
733 'GYFHZIF1,H1Y1XGN,MG', None],
734 ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',
735 'GYFHZI,H1Y1XGN,MG', None],
736 ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z',
737 'GYFLI,I1ZHF1,H1Y1XGN,MG', None]])
738class MCLC(BravaisLattice):
739 conventional_cls = 'MCL'
740 conventional_cellmap = np.array([[1, -1, 0], [1, 1, 0], [0, 0, 1]])
742 def __init__(self, a, b, c, alpha):
743 check_mcl(a, b, c, alpha)
744 super().__init__(a=a, b=b, c=c, alpha=alpha)
746 def _cell(self, a, b, c, alpha):
747 alpha *= np.pi / 180
748 return np.array([[0.5 * a, 0.5 * b, 0], [-0.5 * a, 0.5 * b, 0],
749 [0, c * np.cos(alpha), c * np.sin(alpha)]])
751 def _variant_name(self, a, b, c, alpha):
752 # from ase.geometry.cell import mclc
753 # okay, this is a bit hacky
755 # We need the same parameters here as when determining the points.
756 # Right now we just repeat the code:
757 check_mcl(a, b, c, alpha)
759 a2 = a * a
760 b2 = b * b
761 cosa = np.cos(alpha * _degrees)
762 sina = np.sin(alpha * _degrees)
763 sina2 = sina**2
765 cell = self.tocell()
766 lengths_angles = Cell(cell.reciprocal()).cellpar()
768 kgamma = lengths_angles[-1]
770 eps = self._eps
771 # We should not compare angles in degrees versus lengths with
772 # the same precision.
773 if abs(kgamma - 90) < eps:
774 variant = 2
775 elif kgamma > 90:
776 variant = 1
777 elif kgamma < 90:
778 num = b * cosa / c + b2 * sina2 / a2
779 if abs(num - 1) < eps:
780 variant = 4
781 elif num < 1:
782 variant = 3
783 else:
784 variant = 5
785 variant = 'MCLC' + str(variant)
786 return variant
788 def _special_points(self, a, b, c, alpha, variant):
789 variant = int(variant.name[-1])
791 a2 = a * a
792 b2 = b * b
793 # c2 = c * c
794 cosa = np.cos(alpha * _degrees)
795 sina = np.sin(alpha * _degrees)
796 sina2 = sina**2
798 if variant == 1 or variant == 2:
799 zeta = (2 - b * cosa / c) / (4 * sina2)
800 eta = 0.5 + 2 * zeta * c * cosa / b
801 psi = .75 - a2 / (4 * b2 * sina * sina)
802 phi = psi + (.75 - psi) * b * cosa / c
804 points = [[0, 0, 0],
805 [.5, 0, 0],
806 [0, -.5, 0],
807 [1 - zeta, 1 - zeta, 1 - eta],
808 [zeta, zeta, eta],
809 [-zeta, -zeta, 1 - eta],
810 [1 - zeta, -zeta, 1 - eta],
811 [phi, 1 - phi, .5],
812 [1 - phi, phi - 1, .5],
813 [.5, .5, .5],
814 [.5, 0, .5],
815 [1 - psi, psi - 1, 0],
816 [psi, 1 - psi, 0],
817 [psi - 1, -psi, 0],
818 [.5, .5, 0],
819 [-.5, -.5, 0],
820 [0, 0, .5]]
821 elif variant == 3 or variant == 4:
822 mu = .25 * (1 + b2 / a2)
823 delta = b * c * cosa / (2 * a2)
824 zeta = mu - 0.25 + (1 - b * cosa / c) / (4 * sina2)
825 eta = 0.5 + 2 * zeta * c * cosa / b
826 phi = 1 + zeta - 2 * mu
827 psi = eta - 2 * delta
829 points = [[0, 0, 0],
830 [1 - phi, 1 - phi, 1 - psi],
831 [phi, phi - 1, psi],
832 [1 - phi, -phi, 1 - psi],
833 [zeta, zeta, eta],
834 [1 - zeta, -zeta, 1 - eta],
835 [-zeta, -zeta, 1 - eta],
836 [.5, -.5, .5],
837 [.5, 0, .5],
838 [.5, 0, 0],
839 [0, -.5, 0],
840 [.5, -.5, 0],
841 [mu, mu, delta],
842 [1 - mu, -mu, -delta],
843 [-mu, -mu, -delta],
844 [mu, mu - 1, delta],
845 [0, 0, .5]]
846 elif variant == 5:
847 zeta = .25 * (b2 / a2 + (1 - b * cosa / c) / sina2)
848 eta = 0.5 + 2 * zeta * c * cosa / b
849 mu = .5 * eta + b2 / (4 * a2) - b * c * cosa / (2 * a2)
850 nu = 2 * mu - zeta
851 omega = (4 * nu - 1 - b2 * sina2 / a2) * c / (2 * b * cosa)
852 delta = zeta * c * cosa / b + omega / 2 - .25
853 rho = 1 - zeta * a2 / b2
855 points = [[0, 0, 0],
856 [nu, nu, omega],
857 [1 - nu, 1 - nu, 1 - omega],
858 [nu, nu - 1, omega],
859 [zeta, zeta, eta],
860 [1 - zeta, -zeta, 1 - eta],
861 [-zeta, -zeta, 1 - eta],
862 [rho, 1 - rho, .5],
863 [1 - rho, rho - 1, .5],
864 [.5, .5, .5],
865 [.5, 0, .5],
866 [.5, 0, 0],
867 [0, -.5, 0],
868 [.5, -.5, 0],
869 [mu, mu, delta],
870 [1 - mu, -mu, -delta],
871 [-mu, -mu, -delta],
872 [mu, mu - 1, delta],
873 [0, 0, .5]]
875 return points
878tri_angles_explanation = """\
879Angles kalpha, kbeta and kgamma of TRI lattice must be 1) all greater \
880than 90 degrees with kgamma being the smallest, or 2) all smaller than \
88190 with kgamma being the largest, or 3) kgamma=90 being the \
882smallest of the three, or 4) kgamma=90 being the largest of the three. \
883Angles of reciprocal lattice are kalpha={}, kbeta={}, kgamma={}. \
884If you don't care, please use Cell.fromcellpar() instead."""
886# XXX labels, paths, are all the same.
889@bravaisclass('primitive triclinic', 'triclinic', 'triclinic', 'aP',
890 ('a', 'b', 'c', 'alpha', 'beta', 'gamma'),
891 [['TRI1a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
892 ['TRI2a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
893 ['TRI1b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
894 ['TRI2b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None]])
895class TRI(BravaisLattice):
896 conventional_cls = 'TRI'
897 conventional_cellmap = _identity
899 def __init__(self, a, b, c, alpha, beta, gamma):
900 super().__init__(a=a, b=b, c=c, alpha=alpha, beta=beta,
901 gamma=gamma)
903 def _cell(self, a, b, c, alpha, beta, gamma):
904 alpha, beta, gamma = np.array([alpha, beta, gamma])
905 singamma = np.sin(gamma * _degrees)
906 cosgamma = np.cos(gamma * _degrees)
907 cosbeta = np.cos(beta * _degrees)
908 cosalpha = np.cos(alpha * _degrees)
909 a3x = c * cosbeta
910 a3y = c / singamma * (cosalpha - cosbeta * cosgamma)
911 a3z = c / singamma * np.sqrt(singamma**2 - cosalpha**2 - cosbeta**2
912 + 2 * cosalpha * cosbeta * cosgamma)
913 return np.array([[a, 0, 0], [b * cosgamma, b * singamma, 0],
914 [a3x, a3y, a3z]])
916 def _variant_name(self, a, b, c, alpha, beta, gamma):
917 cell = Cell.new([a, b, c, alpha, beta, gamma])
918 icellpar = Cell(cell.reciprocal()).cellpar()
919 kangles = kalpha, kbeta, kgamma = icellpar[3:]
921 def raise_unconventional():
922 raise UnconventionalLattice(tri_angles_explanation
923 .format(*kangles))
925 eps = self._eps
926 if abs(kgamma - 90) < eps:
927 if kalpha > 90 and kbeta > 90:
928 var = '2a'
929 elif kalpha < 90 and kbeta < 90:
930 var = '2b'
931 else:
932 # Is this possible? Maybe due to epsilon
933 raise_unconventional()
934 elif all(kangles > 90):
935 if kgamma > min(kangles):
936 raise_unconventional()
937 var = '1a'
938 elif all(kangles < 90): # and kgamma > max(kalpha, kbeta):
939 if kgamma < max(kangles):
940 raise_unconventional()
941 var = '1b'
942 else:
943 raise_unconventional()
945 return 'TRI' + var
947 def _special_points(self, a, b, c, alpha, beta, gamma, variant):
948 # (None of the points actually depend on any parameters)
949 # (We should store the points openly on the variant objects)
950 if variant.name == 'TRI1a' or variant.name == 'TRI2a':
951 points = [[0., 0., 0.],
952 [.5, .5, 0],
953 [0, .5, .5],
954 [.5, 0, .5],
955 [.5, .5, .5],
956 [.5, 0, 0],
957 [0, .5, 0],
958 [0, 0, .5]]
959 else:
960 points = [[0, 0, 0],
961 [.5, -.5, 0],
962 [0, 0, .5],
963 [-.5, -.5, .5],
964 [0, -.5, .5],
965 [0, -0.5, 0],
966 [.5, 0, 0],
967 [-.5, 0, .5]]
969 return points
972def get_subset_points(names, points):
973 newpoints = {}
974 for name in names:
975 newpoints[name] = points[name]
977 return newpoints
980@bravaisclass('primitive oblique', 'monoclinic', None, 'mp',
981 ('a', 'b', 'alpha'), [['OBL', 'GYHCH1X', 'GYHCH1XG', None]],
982 ndim=2)
983class OBL(BravaisLattice):
984 def __init__(self, a, b, alpha, **kwargs):
985 check_rect(a, b)
986 if alpha >= 90:
987 raise UnconventionalLattice(
988 f'Expected alpha < 90, got alpha={alpha}')
989 super().__init__(a=a, b=b, alpha=alpha, **kwargs)
991 def _cell(self, a, b, alpha):
992 cosa = np.cos(alpha * _degrees)
993 sina = np.sin(alpha * _degrees)
995 return np.array([[a, 0, 0],
996 [b * cosa, b * sina, 0],
997 [0., 0., 0.]])
999 def _special_points(self, a, b, alpha, variant):
1000 cosa = np.cos(alpha * _degrees)
1001 eta = (1 - a * cosa / b) / (2 * np.sin(alpha * _degrees)**2)
1002 nu = .5 - eta * b * cosa / a
1004 points = [[0, 0, 0],
1005 [0, 0.5, 0],
1006 [eta, 1 - nu, 0],
1007 [.5, .5, 0],
1008 [1 - eta, nu, 0],
1009 [.5, 0, 0]]
1011 return points
1014@bravaisclass('primitive hexagonal', 'hexagonal', None, 'hp', 'a',
1015 [['HEX2D', 'GMK', 'GMKG',
1016 get_subset_points('GMK',
1017 sc_special_points['hexagonal'])]],
1018 ndim=2)
1019class HEX2D(BravaisLattice):
1020 def __init__(self, a, **kwargs):
1021 super().__init__(a=a, **kwargs)
1023 def _cell(self, a):
1024 x = 0.5 * np.sqrt(3)
1025 return np.array([[a, 0, 0],
1026 [-0.5 * a, x * a, 0],
1027 [0., 0., 0.]])
1030def check_rect(a, b):
1031 if a >= b:
1032 raise UnconventionalLattice(f'Expected a < b, got a={a}, b={b}')
1035@bravaisclass('primitive rectangular', 'orthorhombic', None, 'op', 'ab',
1036 [['RECT', 'GXSY', 'GXSYGS',
1037 get_subset_points('GXSY',
1038 sc_special_points['orthorhombic'])]],
1039 ndim=2)
1040class RECT(BravaisLattice):
1041 def __init__(self, a, b, **kwargs):
1042 check_rect(a, b)
1043 super().__init__(a=a, b=b, **kwargs)
1045 def _cell(self, a, b):
1046 return np.array([[a, 0, 0],
1047 [0, b, 0],
1048 [0, 0, 0.]])
1051@bravaisclass('centred rectangular', 'orthorhombic', None, 'oc',
1052 ('a', 'alpha'), [['CRECT', 'GXA1Y', 'GXA1YG', None]], ndim=2)
1053class CRECT(BravaisLattice):
1054 def __init__(self, a, alpha, **kwargs):
1055 # It would probably be better to define the CRECT cell
1056 # by (a, b) rather than (a, alpha). Then we can require a < b
1057 # like in ordinary RECT.
1058 #
1059 # In 3D, all lattices in the same family generally take
1060 # identical parameters.
1061 if alpha >= 90:
1062 raise UnconventionalLattice(
1063 f'Expected alpha < 90. Got alpha={alpha}')
1064 super().__init__(a=a, alpha=alpha, **kwargs)
1066 def _cell(self, a, alpha):
1067 x = np.cos(alpha * _degrees)
1068 y = np.sin(alpha * _degrees)
1069 return np.array([[a, 0, 0],
1070 [a * x, a * y, 0],
1071 [0, 0, 0.]])
1073 def _special_points(self, a, alpha, variant):
1074 sina2 = np.sin(alpha / 2 * _degrees)**2
1075 sina = np.sin(alpha * _degrees)**2
1076 eta = sina2 / sina
1077 cosa = np.cos(alpha * _degrees)
1078 xi = eta * cosa
1080 points = [[0, 0, 0],
1081 [eta, - eta, 0],
1082 [0.5 + xi, 0.5 - xi, 0],
1083 [0.5, 0.5, 0]]
1085 return points
1088@bravaisclass('primitive square', 'tetragonal', None, 'tp', ('a',),
1089 [['SQR', 'GMX', 'MGXM',
1090 get_subset_points('GMX', sc_special_points['tetragonal'])]],
1091 ndim=2)
1092class SQR(BravaisLattice):
1093 def __init__(self, a, **kwargs):
1094 super().__init__(a=a, **kwargs)
1096 def _cell(self, a):
1097 return np.array([[a, 0, 0],
1098 [0, a, 0],
1099 [0, 0, 0.]])
1102@bravaisclass('primitive line', 'line', None, '?', ('a',),
1103 [['LINE', 'GX', 'GX', {'G': [0, 0, 0], 'X': [0.5, 0, 0]}]],
1104 ndim=1)
1105class LINE(BravaisLattice):
1106 def __init__(self, a, **kwargs):
1107 super().__init__(a=a, **kwargs)
1109 def _cell(self, a):
1110 return np.array([[a, 0.0, 0.0],
1111 [0.0, 0.0, 0.0],
1112 [0.0, 0.0, 0.0]])
1115def celldiff(cell1, cell2):
1116 """Return a unitless measure of the difference between two cells."""
1117 cell1 = Cell.ascell(cell1).complete()
1118 cell2 = Cell.ascell(cell2).complete()
1119 v1v2 = cell1.volume * cell2.volume
1120 if v1v2 < 1e-10:
1121 # (Proposed cell may be linearly dependent)
1122 return np.inf
1124 scale = v1v2**(-1. / 3.) # --> 1/Ang^2
1125 x1 = cell1 @ cell1.T
1126 x2 = cell2 @ cell2.T
1127 dev = scale * np.abs(x2 - x1).max()
1128 return dev
1131def get_lattice_from_canonical_cell(cell, eps=2e-4):
1132 """Return a Bravais lattice representing the given cell.
1134 This works only for cells that are derived from the standard form
1135 (as generated by lat.tocell()) or rotations thereof.
1137 If the given cell does not resemble the known form of a Bravais
1138 lattice, raise RuntimeError."""
1139 return LatticeChecker(cell, eps).match()
1142def identify_lattice(cell, eps=2e-4, *, pbc=True):
1143 """Find Bravais lattice representing this cell.
1145 Returns Bravais lattice object representing the cell along with
1146 an operation that, applied to the cell, yields the same lengths
1147 and angles as the Bravais lattice object."""
1148 from ase.geometry.bravais_type_engine import niggli_op_table
1150 pbc = cell.any(1) & pbc2pbc(pbc)
1151 npbc = sum(pbc)
1153 cell = cell.uncomplete(pbc)
1154 rcell, reduction_op = cell.niggli_reduce(eps=eps)
1156 # We tabulate the cell's Niggli-mapped versions so we don't need to
1157 # redo any work when the same Niggli-operation appears multiple times
1158 # in the table:
1159 memory = {}
1161 # We loop through the most symmetric kinds (CUB etc.) and return
1162 # the first one we find:
1163 for latname in LatticeChecker.check_orders[npbc]:
1164 # There may be multiple Niggli operations that produce valid
1165 # lattices, at least for MCL. In that case we will pick the
1166 # one whose angle is closest to 90, but it means we cannot
1167 # just return the first one we find so we must remember then:
1168 matching_lattices = []
1170 for op_key in niggli_op_table[latname]:
1171 checker_and_op = memory.get(op_key)
1172 if checker_and_op is None:
1173 normalization_op = np.array(op_key).reshape(3, 3)
1174 candidate = Cell(np.linalg.inv(normalization_op.T) @ rcell)
1175 checker = LatticeChecker(candidate, eps=eps)
1176 memory[op_key] = (checker, normalization_op)
1177 else:
1178 checker, normalization_op = checker_and_op
1180 lat = checker.query(latname)
1181 if lat is not None:
1182 op = normalization_op @ np.linalg.inv(reduction_op)
1183 matching_lattices.append((lat, op))
1185 # Among any matching lattices, return the one with lowest
1186 # orthogonality defect:
1187 best = None
1188 best_defect = np.inf
1189 for lat, op in matching_lattices:
1190 cell = lat.tocell()
1191 lengths = cell.lengths()[pbc]
1192 generalized_volume = cell.complete().volume
1193 defect = np.prod(lengths) / generalized_volume
1194 if defect < best_defect:
1195 best = lat, op
1196 best_defect = defect
1198 if best is not None:
1199 return best
1201 raise RuntimeError('Failed to recognize lattice')
1204class LatticeChecker:
1205 # The check order is slightly different than elsewhere listed order
1206 # as we need to check HEX/RHL before the ORCx family.
1207 check_orders = {
1208 1: ['LINE'],
1209 2: ['SQR', 'RECT', 'HEX2D', 'CRECT', 'OBL'],
1210 3: ['CUB', 'FCC', 'BCC', 'TET', 'BCT', 'HEX', 'RHL',
1211 'ORC', 'ORCF', 'ORCI', 'ORCC', 'MCL', 'MCLC', 'TRI']}
1213 def __init__(self, cell, eps=2e-4):
1214 """Generate Bravais lattices that look (or not) like the given cell.
1216 The cell must be reduced to canonical form, i.e., it must
1217 be possible to produce a cell with the same lengths and angles
1218 by directly through one of the Bravais lattice classes.
1220 Generally for internal use (this module).
1222 For each of the 14 Bravais lattices, this object can produce
1223 a lattice object which represents the same cell, or None if
1224 the tolerance eps is not met."""
1225 self.cell = cell
1226 self.eps = eps
1228 self.cellpar = cell.cellpar()
1229 self.lengths = self.A, self.B, self.C = self.cellpar[:3]
1230 self.angles = self.cellpar[3:]
1232 # Use a 'neutral' length for checking cubic lattices
1233 self.A0 = self.lengths.mean()
1235 # Vector of the diagonal and then off-diagonal dot products:
1236 # [a1 · a1, a2 · a2, a3 · a3, a2 · a3, a3 · a1, a1 · a2]
1237 self.prods = (cell @ cell.T).flat[[0, 4, 8, 5, 2, 1]]
1239 def _check(self, latcls, *args):
1240 if any(arg <= 0 for arg in args):
1241 return None
1242 try:
1243 lat = latcls(*args)
1244 except UnconventionalLattice:
1245 return None
1247 newcell = lat.tocell()
1248 err = celldiff(self.cell, newcell)
1249 if err < self.eps:
1250 return lat
1252 def match(self):
1253 """Match cell against all lattices, returning most symmetric match.
1255 Returns the lattice object. Raises RuntimeError on failure."""
1256 for name in self.check_orders[self.cell.rank]:
1257 lat = self.query(name)
1258 if lat:
1259 return lat
1260 else:
1261 raise RuntimeError('Could not find lattice type for cell '
1262 'with lengths and angles {}'
1263 .format(self.cell.cellpar().tolist()))
1265 def query(self, latname):
1266 """Match cell against named Bravais lattice.
1268 Return lattice object on success, None on failure."""
1269 meth = getattr(self, latname)
1270 lat = meth()
1271 return lat
1273 def LINE(self):
1274 return self._check(LINE, self.lengths[0])
1276 def SQR(self):
1277 return self._check(SQR, self.lengths[0])
1279 def RECT(self):
1280 return self._check(RECT, *self.lengths[:2])
1282 def CRECT(self):
1283 return self._check(CRECT, self.lengths[0], self.angles[2])
1285 def HEX2D(self):
1286 return self._check(HEX2D, self.lengths[0])
1288 def OBL(self):
1289 return self._check(OBL, *self.lengths[:2], self.angles[2])
1291 def CUB(self):
1292 # These methods (CUB, FCC, ...) all return a lattice object if
1293 # it matches, else None.
1294 return self._check(CUB, self.A0)
1296 def FCC(self):
1297 return self._check(FCC, np.sqrt(2) * self.A0)
1299 def BCC(self):
1300 return self._check(BCC, 2.0 * self.A0 / np.sqrt(3))
1302 def TET(self):
1303 return self._check(TET, self.A, self.C)
1305 def _bct_orci_lengths(self):
1306 # Coordinate-system independent relation for BCT and ORCI
1307 # standard cells:
1308 # a1 · a1 + a2 · a3 == a² / 2
1309 # a2 · a2 + a3 · a1 == a² / 2 (BCT)
1310 # == b² / 2 (ORCI)
1311 # a3 · a3 + a1 · a2 == c² / 2
1312 # We use these to get a, b, and c in those cases.
1313 prods = self.prods
1314 lengthsqr = 2.0 * (prods[:3] + prods[3:])
1315 if any(lengthsqr < 0):
1316 return None
1317 return np.sqrt(lengthsqr)
1319 def BCT(self):
1320 lengths = self._bct_orci_lengths()
1321 if lengths is None:
1322 return None
1323 return self._check(BCT, lengths[0], lengths[2])
1325 def HEX(self):
1326 return self._check(HEX, self.A, self.C)
1328 def RHL(self):
1329 return self._check(RHL, self.A, self.angles[0])
1331 def ORC(self):
1332 return self._check(ORC, *self.lengths)
1334 def ORCF(self):
1335 # ORCF standard cell:
1336 # a2 · a3 = a²/4
1337 # a3 · a1 = b²/4
1338 # a1 · a2 = c²/4
1339 prods = self.prods
1340 if all(prods[3:] > 0):
1341 orcf_abc = 2 * np.sqrt(prods[3:])
1342 return self._check(ORCF, *orcf_abc)
1344 def ORCI(self):
1345 lengths = self._bct_orci_lengths()
1346 if lengths is None:
1347 return None
1348 return self._check(ORCI, *lengths)
1350 def _orcc_ab(self):
1351 # ORCC: a1 · a1 + a2 · a3 = a²/2
1352 # a2 · a2 - a2 · a3 = b²/2
1353 prods = self.prods
1354 orcc_sqr_ab = np.empty(2)
1355 orcc_sqr_ab[0] = 2.0 * (prods[0] + prods[5])
1356 orcc_sqr_ab[1] = 2.0 * (prods[1] - prods[5])
1357 if all(orcc_sqr_ab > 0):
1358 return np.sqrt(orcc_sqr_ab)
1360 def ORCC(self):
1361 orcc_lengths_ab = self._orcc_ab()
1362 if orcc_lengths_ab is None:
1363 return None
1364 return self._check(ORCC, *orcc_lengths_ab, self.C)
1366 def MCL(self):
1367 return self._check(MCL, *self.lengths, self.angles[0])
1369 def MCLC(self):
1370 # MCLC is similar to ORCC:
1371 orcc_ab = self._orcc_ab()
1372 if orcc_ab is None:
1373 return None
1375 prods = self.prods
1376 C = self.C
1377 mclc_a, mclc_b = orcc_ab[::-1] # a, b reversed wrt. ORCC
1378 mclc_cosa = 2.0 * prods[3] / (mclc_b * C)
1379 if -1 < mclc_cosa < 1:
1380 mclc_alpha = np.arccos(mclc_cosa) * 180 / np.pi
1381 if mclc_b > C:
1382 # XXX Temporary fix for certain otherwise
1383 # unrecognizable lattices.
1384 #
1385 # This error could happen if the input lattice maps to
1386 # something just outside the domain of conventional
1387 # lattices (less than the tolerance). Our solution is to
1388 # propose a nearby conventional lattice instead, which
1389 # will then be accepted if it's close enough.
1390 mclc_b = 0.5 * (mclc_b + C)
1391 C = mclc_b
1392 return self._check(MCLC, mclc_a, mclc_b, C, mclc_alpha)
1394 def TRI(self):
1395 return self._check(TRI, *self.cellpar)
1398def all_variants():
1399 """For testing and examples; yield all variants of all lattices."""
1400 a, b, c = 3., 4., 5.
1401 alpha = 55.0
1402 yield CUB(a)
1403 yield FCC(a)
1404 yield BCC(a)
1405 yield TET(a, c)
1406 bct1 = BCT(2 * a, c)
1407 bct2 = BCT(a, c)
1408 assert bct1.variant == 'BCT1'
1409 assert bct2.variant == 'BCT2'
1411 yield bct1
1412 yield bct2
1414 yield ORC(a, b, c)
1416 a0 = np.sqrt(1.0 / (1 / b**2 + 1 / c**2))
1417 orcf1 = ORCF(0.5 * a0, b, c)
1418 orcf2 = ORCF(1.2 * a0, b, c)
1419 orcf3 = ORCF(a0, b, c)
1420 assert orcf1.variant == 'ORCF1'
1421 assert orcf2.variant == 'ORCF2'
1422 assert orcf3.variant == 'ORCF3'
1423 yield orcf1
1424 yield orcf2
1425 yield orcf3
1427 yield ORCI(a, b, c)
1428 yield ORCC(a, b, c)
1430 yield HEX(a, c)
1432 rhl1 = RHL(a, alpha=55.0)
1433 assert rhl1.variant == 'RHL1'
1434 yield rhl1
1436 rhl2 = RHL(a, alpha=105.0)
1437 assert rhl2.variant == 'RHL2'
1438 yield rhl2
1440 # With these lengths, alpha < 65 (or so) would result in a lattice that
1441 # could also be represented with alpha > 65, which is more conventional.
1442 yield MCL(a, b, c, alpha=70.0)
1444 mclc1 = MCLC(a, b, c, 80)
1445 assert mclc1.variant == 'MCLC1'
1446 yield mclc1
1447 # mclc2 has same special points as mclc1
1449 mclc3 = MCLC(1.8 * a, b, c * 2, 80)
1450 assert mclc3.variant == 'MCLC3'
1451 yield mclc3
1452 # mclc4 has same special points as mclc3
1454 # XXX We should add MCLC2 and MCLC4 as well.
1456 mclc5 = MCLC(b, b, 1.1 * b, 70)
1457 assert mclc5.variant == 'MCLC5'
1458 yield mclc5
1460 def get_tri(kcellpar):
1461 # We build the TRI lattices from cellpars of reciprocal cell
1462 icell = Cell.fromcellpar(kcellpar)
1463 cellpar = Cell(4 * icell.reciprocal()).cellpar()
1464 return TRI(*cellpar)
1466 tri1a = get_tri([1., 1.2, 1.4, 120., 110., 100.])
1467 assert tri1a.variant == 'TRI1a'
1468 yield tri1a
1470 tri1b = get_tri([1., 1.2, 1.4, 50., 60., 70.])
1471 assert tri1b.variant == 'TRI1b'
1472 yield tri1b
1474 tri2a = get_tri([1., 1.2, 1.4, 120., 110., 90.])
1475 assert tri2a.variant == 'TRI2a'
1476 yield tri2a
1477 tri2b = get_tri([1., 1.2, 1.4, 50., 60., 90.])
1478 assert tri2b.variant == 'TRI2b'
1479 yield tri2b
1481 # Choose an OBL lattice that round-trip-converts to itself.
1482 # The default a/b/alpha parameters result in another representation
1483 # of the same lattice.
1484 yield OBL(a=3.0, b=3.35, alpha=77.85)
1485 yield RECT(a, b)
1486 yield CRECT(a, alpha=alpha)
1487 yield HEX2D(a)
1488 yield SQR(a)
1489 yield LINE(a)