1from abc import abstractmethod, ABC
2import functools
3import warnings
4import numpy as np
5from typing import Dict, List
6
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
11
12
13@functools.wraps(newbulk)
14def bulk(*args, **kwargs):
15    warnings.warn('Use ase.build.bulk() instead', stacklevel=2)
16    return newbulk(*args, **kwargs)
17
18
19_degrees = np.pi / 180
20
21
22class BravaisLattice(ABC):
23    """Represent Bravais lattices and data related to the Brillouin zone.
24
25    There are 14 3D Bravais classes: CUB, FCC, BCC, ..., and TRI, and
26    five 2D classes.
27
28    Each class stores basic static information:
29
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')
39
40    Each class can be instantiated with the specific lattice parameters
41    that apply to that lattice:
42
43    >>> MCL(3, 4, 5, 80)
44    MCL(a=3, b=4, c=5, alpha=80)
45
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
55
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
64
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]
71
72    @property
73    def variant(self) -> str:
74        """Return name of lattice variant.
75
76        >>> BCT(3, 5).variant
77        'BCT2'
78        """
79        return self._variant.name
80
81    def __getattr__(self, name: str):
82        if name in self._parameters:
83            return self._parameters[name]
84        return self.__getattribute__(name)  # Raises error
85
86    def vars(self) -> Dict[str, float]:
87        """Get parameter names and values of this lattice as a dictionary."""
88        return dict(self._parameters)
89
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)
94
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)
99
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
108
109    def cellpar(self) -> np.ndarray:
110        """Get cell lengths and angles as array of length 6.
111
112        See :func:`ase.geometry.Cell.cellpar`."""
113        # (Just a brute-force implementation)
114        cell = self.tocell()
115        return cell.cellpar()
116
117    @property
118    def special_path(self) -> str:
119        """Get default special k-point path for this lattice as a string.
120
121        >>> BCT(3, 5).special_path
122        'GXYSGZS1NPY1Z,XP'
123        """
124        return self._variant.special_path
125
126    @property
127    def special_point_names(self) -> List[str]:
128        """Return all special point names as a list of strings.
129
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]
136
137    def get_special_points_array(self) -> np.ndarray:
138        """Return all special points for this lattice as an array.
139
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
150
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)
156
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
161
162        labels = self.special_point_names
163        points = self.get_special_points_array()
164
165        return dict(zip(labels, points))
166
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)
173
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.
177
178        See :meth:`ase.cell.Cell.bandpath` for description of parameters.
179
180        >>> BCT(3, 5).bandpath()
181        BandPath(path='GXYSGZS1NPY1Z,XP', cell=[3x3], \
182special_points={GNPSS1XYY1Z}, kpts=[51x3])
183
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.
188
189        """
190        if special_points is None:
191            special_points = self.get_special_points()
192
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)
200
201        cell = self.tocell()
202        if transformation is not None:
203            cell = transformation.dot(cell)
204
205        bandpath = BandPath(cell=cell, path=path,
206                            special_points=special_points)
207        return bandpath.interpolate(npoints=npoints, density=density)
208
209    @abstractmethod
210    def _cell(self, **kwargs):
211        """Return a Cell object from this Bravais lattice.
212
213        Arguments are the dictionary of Bravais parameters."""
214        pass
215
216    def _special_points(self, **kwargs):
217        """Return the special point coordinates as an npoints x 3 sequence.
218
219        Subclasses typically return a nested list.
220
221        Ordering must be same as kpoint labels.
222
223        Arguments are the dictionary of Bravais parameters and the variant."""
224        raise NotImplementedError
225
226    def _variant_name(self, **kwargs):
227        """Return the name (e.g. 'ORCF3') of variant.
228
229        Arguments will be the dictionary of Bravais parameters."""
230        raise NotImplementedError
231
232    def __format__(self, spec):
233        tokens = []
234        if not spec:
235            spec = '.6g'
236        template = '{}={:%s}' % spec
237
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))
242
243    def __str__(self) -> str:
244        return self.__format__('')
245
246    def __repr__(self) -> str:
247        return self.__format__('.20g')
248
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
253
254        coordstring = '\n'.join(['    {:2s} {:7.4f} {:7.4f} {:7.4f}'
255                                 .format(label, *points[label])
256                                 for label in labels])
257
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
267
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))
276
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)
284
285        return '\n'.join(chunks)
286
287
288class Variant:
289    variant_desc = """\
290Variant name: {name}
291  Special point names: {special_point_names}
292  Default path: {special_path}
293"""
294
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)
307
308    def __str__(self) -> str:
309        return self.variant_desc.format(**vars(self))
310
311
312bravais_names = []
313bravais_lattices = {}
314bravais_classes = {}
315
316
317def bravaisclass(longname, crystal_family, lattice_system, pearson_symbol,
318                 parameters, variants, ndim=3):
319    """Decorator for Bravais lattice classes.
320
321    This sets a number of class variables and processes the information
322    about different variants into a list of Variant objects."""
323
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
335
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)
341
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
347
348    return decorate
349
350
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]])
355
356
357class UnconventionalLattice(ValueError):
358    pass
359
360
361class Cubic(BravaisLattice):
362    """Abstract class for cubic lattices."""
363    conventional_cls = 'CUB'
364
365    def __init__(self, a):
366        BravaisLattice.__init__(self, a=a)
367
368
369@bravaisclass('primitive cubic', 'cubic', 'cubic', 'cP', 'a',
370              [['CUB', 'GXRM', 'GXMGRX,MR', sc_special_points['cubic']]])
371class CUB(Cubic):
372    conventional_cellmap = _identity
373
374    def _cell(self, a):
375        return a * np.eye(3)
376
377
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
382
383    def _cell(self, a):
384        return 0.5 * np.array([[0., a, a], [a, 0, a], [a, a, 0]])
385
386
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
391
392    def _cell(self, a):
393        return 0.5 * np.array([[-a, a, a], [a, -a, a], [a, a, -a]])
394
395
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
402
403    def __init__(self, a, c):
404        BravaisLattice.__init__(self, a=a, c=c)
405
406    def _cell(self, a, c):
407        return np.diag(np.array([a, a, c]))
408
409
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
419
420    def __init__(self, a, c):
421        BravaisLattice.__init__(self, a=a, c=c)
422
423    def _cell(self, a, c):
424        return 0.5 * np.array([[-a, a, c], [a, -a, c], [a, a, -c]])
425
426    def _variant_name(self, a, c):
427        return 'BCT1' if c < a else 'BCT2'
428
429    def _special_points(self, a, c, variant):
430        a2 = a * a
431        c2 = c * c
432
433        assert variant.name in self.variants
434
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
457
458
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))
463
464
465class Orthorhombic(BravaisLattice):
466    """Abstract class for orthorhombic types."""
467
468    def __init__(self, a, b, c):
469        check_orc(a, b, c)
470        BravaisLattice.__init__(self, a=a, b=b, c=c)
471
472
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
480
481    def _cell(self, a, b, c):
482        return np.diag([a, b, c]).astype(float)
483
484
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
493
494    def _cell(self, a, b, c):
495        return 0.5 * np.array([[0, b, c], [a, 0, c], [a, b, 0]])
496
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)
503
504        if variant.name == 'ORCF1' or variant.name == 'ORCF3':
505            zeta = xminus
506            eta = xplus
507
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
522
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.]]
534
535        return points
536
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'
542
543
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
550
551    def _cell(self, a, b, c):
552        return 0.5 * np.array([[-a, b, c], [a, -b, c], [a, b, -c]])
553
554    def _special_points(self, a, b, c, variant):
555        a2 = a**2
556        b2 = b**2
557        c2 = c**2
558
559        zeta = .25 * (1 + a2 / c2)
560        eta = .25 * (1 + b2 / c2)
561        delta = .25 * (b2 - a2) / c2
562        mu = .25 * (a2 + b2) / c2
563
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
578
579
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]])
586
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)
592
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]])
596
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
610
611
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
619
620    def __init__(self, a, c):
621        BravaisLattice.__init__(self, a=a, c=c)
622
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]])
627
628
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
636
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)
642
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]])
653
654    def _variant_name(self, a, alpha):
655        return 'RHL1' if alpha < 90 else 'RHL2'
656
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
686
687
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))
693
694
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
701
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)
705
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)]])
710
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
715
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
733
734    def _variant_name(self, a, b, c, alpha):
735        check_mcl(a, b, c, alpha)
736        return 'MCL'
737
738
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]])
754
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)
758
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)]])
763
764    def _variant_name(self, a, b, c, alpha):
765        # from ase.geometry.cell import mclc
766        # okay, this is a bit hacky
767
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)
771
772        a2 = a * a
773        b2 = b * b
774        cosa = np.cos(alpha * _degrees)
775        sina = np.sin(alpha * _degrees)
776        sina2 = sina**2
777
778        cell = self.tocell()
779        lengths_angles = Cell(cell.reciprocal()).cellpar()
780
781        kgamma = lengths_angles[-1]
782
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
800
801    def _special_points(self, a, b, c, alpha, variant):
802        variant = int(variant.name[-1])
803
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
810
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
816
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
841
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
867
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]]
887
888        return points
889
890
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."""
898
899# XXX labels, paths, are all the same.
900
901
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
911
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)
915
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]])
928
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:]
933
934        def raise_unconventional():
935            raise UnconventionalLattice(tri_angles_explanation
936                                        .format(*kangles))
937
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()
957
958        return 'TRI' + var
959
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]]
981
982        return points
983
984
985def get_subset_points(names, points):
986    newpoints = {}
987    for name in names:
988        newpoints[name] = points[name]
989
990    return newpoints
991
992
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)
999
1000    def _cell(self, a, b, alpha):
1001        cosa = np.cos(alpha * _degrees)
1002        sina = np.sin(alpha * _degrees)
1003
1004        return np.array([[a, 0, 0],
1005                         [b * cosa, b * sina, 0],
1006                         [0., 0., 0.]])
1007
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
1015
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
1019
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]]
1026
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)
1033
1034        return points
1035
1036
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)
1045
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.]])
1051
1052
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)
1061
1062    def _cell(self, a, b):
1063        return np.array([[a, 0, 0],
1064                         [0, b, 0],
1065                         [0, 0, 0.]])
1066
1067
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)
1073
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.]])
1080
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
1091
1092        points = [[0, 0, 0],
1093                  [eta, - eta, 0],
1094                  [0.5 + xi, 0.5 - xi, 0],
1095                  [0.5, 0.5, 0]]
1096
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
1104
1105
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)
1113
1114    def _cell(self, a):
1115        return np.array([[a, 0, 0],
1116                         [0, a, 0],
1117                         [0, 0, 0.]])
1118
1119
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)
1126
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]])
1131
1132
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
1145
1146
1147def get_lattice_from_canonical_cell(cell, eps=2e-4):
1148    """Return a Bravais lattice representing the given cell.
1149
1150    This works only for cells that are derived from the standard form
1151    (as generated by lat.tocell()) or rotations thereof.
1152
1153    If the given cell does not resemble the known form of a Bravais
1154    lattice, raise RuntimeError."""
1155    return LatticeChecker(cell, eps).match()
1156
1157
1158def identify_lattice(cell, eps=2e-4, *, pbc=True):
1159    """Find Bravais lattice representing this cell.
1160
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."""
1164
1165    pbc = cell.any(1) & pbc2pbc(pbc)
1166    npbc = sum(pbc)
1167
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
1181
1182    if npbc == 2:
1183        lat, op = get_2d_bravais_lattice(cell, eps, pbc=pbc)
1184        return lat, op
1185
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))
1192
1193    from ase.geometry.bravais_type_engine import niggli_op_table
1194
1195    if cell.rank < 3:
1196        raise ValueError('Expected 3 linearly independent cell vectors')
1197    rcell, reduction_op = cell.niggli_reduce(eps=eps)
1198
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 = {}
1203
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 = []
1212
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
1222
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))
1227
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
1239
1240        if best is not None:
1241            return best
1242
1243
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']
1249
1250    def __init__(self, cell, eps=2e-4):
1251        """Generate Bravais lattices that look (or not) like the given cell.
1252
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.
1256
1257        Generally for internal use (this module).
1258
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
1264
1265        self.cellpar = cell.cellpar()
1266        self.lengths = self.A, self.B, self.C = self.cellpar[:3]
1267        self.angles = self.cellpar[3:]
1268
1269        # Use a 'neutral' length for checking cubic lattices
1270        self.A0 = self.lengths.mean()
1271
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]]
1275
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
1283
1284        newcell = lat.tocell()
1285        err = celldiff(self.cell, newcell)
1286        if err < self.eps:
1287            return lat
1288
1289    def match(self):
1290        """Match cell against all lattices, returning most symmetric match.
1291
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()))
1301
1302    def query(self, latname):
1303        """Match cell against named Bravais lattice.
1304
1305        Return lattice object on success, None on failure."""
1306        meth = getattr(self, latname)
1307        lat = meth()
1308        return lat
1309
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)
1314
1315    def FCC(self):
1316        return self._check(FCC, np.sqrt(2) * self.A0)
1317
1318    def BCC(self):
1319        return self._check(BCC, 2.0 * self.A0 / np.sqrt(3))
1320
1321    def TET(self):
1322        return self._check(TET, self.A, self.C)
1323
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)
1337
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])
1343
1344    def HEX(self):
1345        return self._check(HEX, self.A, self.C)
1346
1347    def RHL(self):
1348        return self._check(RHL, self.A, self.angles[0])
1349
1350    def ORC(self):
1351        return self._check(ORC, *self.lengths)
1352
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)
1362
1363    def ORCI(self):
1364        lengths = self._bct_orci_lengths()
1365        if lengths is None:
1366            return None
1367        return self._check(ORCI, *lengths)
1368
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)
1378
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)
1384
1385    def MCL(self):
1386        return self._check(MCL, *self.lengths, self.angles[0])
1387
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
1393
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)
1412
1413    def TRI(self):
1414        return self._check(TRI, *self.cellpar)
1415
1416
1417class UnsupportedLattice(ValueError):
1418    pass
1419
1420
1421def get_2d_bravais_lattice(origcell, eps=2e-4, *, pbc=True):
1422
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))
1427
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))
1440
1441    def allclose(a, b):
1442        return np.allclose(a, b, atol=eps)
1443
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]
1452
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
1457
1458        all_lengths_equal = abs(a - b) < eps
1459
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
1477
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
1488
1489    return finallat, finalop.T
1490
1491
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'
1504
1505    yield bct1
1506    yield bct2
1507
1508    yield ORC(a, b, c)
1509
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
1520
1521    yield ORCI(a, b, c)
1522    yield ORCC(a, b, c)
1523
1524    yield HEX(a, c)
1525
1526    rhl1 = RHL(a, alpha=55.0)
1527    assert rhl1.variant == 'RHL1'
1528    yield rhl1
1529
1530    rhl2 = RHL(a, alpha=105.0)
1531    assert rhl2.variant == 'RHL2'
1532    yield rhl2
1533
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)
1537
1538    mclc1 = MCLC(a, b, c, 80)
1539    assert mclc1.variant == 'MCLC1'
1540    yield mclc1
1541    # mclc2 has same special points as mclc1
1542
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
1547
1548    # XXX We should add MCLC2 and MCLC4 as well.
1549
1550    mclc5 = MCLC(b, b, 1.1 * b, 70)
1551    assert mclc5.variant == 'MCLC5'
1552    yield mclc5
1553
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)
1559
1560    tri1a = get_tri([1., 1.2, 1.4, 120., 110., 100.])
1561    assert tri1a.variant == 'TRI1a'
1562    yield tri1a
1563
1564    tri1b = get_tri([1., 1.2, 1.4, 50., 60., 70.])
1565    assert tri1b.variant == 'TRI1b'
1566    yield tri1b
1567
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
1574
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)
1581
1582    if include_blunt_angles:
1583        beta = 110
1584        yield OBL(a, b, alpha=beta)
1585        yield CRECT(a, alpha=beta)
1586