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