Source code for ase.lattice

# flake8: noqa
from abc import abstractmethod, ABC
import functools
import warnings
import numpy as np

from ase.cell import Cell
from ase.build.bulk import bulk as newbulk
from ase.dft.kpoints import parse_path_string, sc_special_points, BandPath
from ase.utils import pbc2pbc


@functools.wraps(newbulk)
def bulk(*args, **kwargs):
    warnings.warn('Use ase.build.bulk() instead', stacklevel=2)
    return newbulk(*args, **kwargs)


_degrees = np.pi / 180


[docs]class BravaisLattice(ABC): """Represent Bravais lattices and data related to the Brillouin zone. There are 14 3D Bravais classes: CUB, FCC, BCC, ..., and TRI, and five 2D classes. Each class stores basic static information: >>> from ase.lattice import FCC, MCL >>> FCC.name 'FCC' >>> FCC.longname 'face-centred cubic' >>> FCC.pearson_symbol 'cF' >>> MCL.parameters ('a', 'b', 'c', 'alpha') Each class can be instantiated with the specific lattice parameters that apply to that lattice: >>> MCL(3, 4, 5, 80) MCL(a=3, b=4, c=5, alpha=80) """ # These parameters can be set by the @bravais decorator for a subclass. # (We could also use metaclasses to do this, but that's more abstract) name = None # e.g. 'CUB', 'BCT', 'ORCF', ... longname = None # e.g. 'cubic', 'body-centred tetragonal', ... parameters = None # e.g. ('a', 'c') variants = None # e.g. {'BCT1': <variant object>, # 'BCT2': <variant object>} ndim = None def __init__(self, **kwargs): p = {} eps = kwargs.pop('eps', 2e-4) for k, v in kwargs.items(): p[k] = float(v) assert set(p) == set(self.parameters) self._parameters = p self._eps = eps if len(self.variants) == 1: # If there's only one it has the same name as the lattice type self._variant = self.variants[self.name] else: name = self._variant_name(**self._parameters) self._variant = self.variants[name] @property def variant(self): """Return name of lattice variant. >>> BCT(3, 5).variant 'BCT2' """ return self._variant.name def __getattr__(self, name): if name in self._parameters: return self._parameters[name] return self.__getattribute__(name) # Raises error
[docs] def vars(self): """Get parameter names and values of this lattice as a dictionary.""" return dict(self._parameters)
[docs] def tocell(self): """Return this lattice as a :class:`~ase.cell.Cell` object.""" cell = self._cell(**self._parameters) return Cell(cell)
def get_transformation(self, cell): # Get transformation matrix relating input cell to canonical cell T = cell.dot(np.linalg.pinv(self.tocell())) msg = 'This transformation changes the length/area/volume of the cell' assert np.isclose(np.abs(np.linalg.det(T[:self.ndim, :self.ndim])), 1), msg return T
[docs] def cellpar(self): """Get cell lengths and angles as array of length 6. See :func:`ase.geometry.Cell.cellpar`.""" # (Just a brute-force implementation) cell = self.tocell() return cell.cellpar()
@property def special_path(self): """Get default special k-point path for this lattice as a string. >>> BCT(3, 5).special_path 'GXYSGZS1NPY1Z,XP' """ return self._variant.special_path @property def special_point_names(self): """Return all special point names as a list of strings. >>> BCT(3, 5).special_point_names ['G', 'N', 'P', 'S', 'S1', 'X', 'Y', 'Y1', 'Z'] """ labels = parse_path_string(self._variant.special_point_names) assert len(labels) == 1 # list of lists return labels[0]
[docs] def get_special_points_array(self): """Return all special points for this lattice as an array. Ordering is consistent with special_point_names.""" if self._variant.special_points is not None: # Fixed dictionary of special points d = self.get_special_points() labels = self.special_point_names assert len(d) == len(labels) points = np.empty((len(d), 3)) for i, label in enumerate(labels): points[i] = d[label] return points # Special points depend on lattice parameters: points = self._special_points(variant=self._variant, **self._parameters) assert len(points) == len(self.special_point_names) return np.array(points)
[docs] def get_special_points(self): """Return a dictionary of named special k-points for this lattice.""" if self._variant.special_points is not None: return self._variant.special_points labels = self.special_point_names points = self.get_special_points_array() return dict(zip(labels, points))
[docs] def plot_bz(self, path=None, special_points=None, **plotkwargs): """Plot the reciprocal cell and default bandpath.""" # Create a generic bandpath (no interpolated kpoints): bandpath = self.bandpath(path=path, special_points=special_points, npoints=0) return bandpath.plot(dimension=self.ndim, **plotkwargs)
[docs] def bandpath(self, path=None, npoints=None, special_points=None, density=None, transformation=None): """Return a :class:`~ase.dft.kpoints.BandPath` for this lattice. See :meth:`ase.cell.Cell.bandpath` for description of parameters. >>> BCT(3, 5).bandpath() BandPath(path='GXYSGZS1NPY1Z,XP', cell=[3x3], special_points={GNPSS1XYY1Z}, kpts=[51x3]) .. note:: This produces the standard band path following AFlow conventions. If your cell does not follow this convention, you will need :meth:`ase.cell.Cell.bandpath` instead or the kpoints may not correspond to your particular cell. """ if special_points is None: special_points = self.get_special_points() if path is None: path = self._variant.special_path elif not isinstance(path, str): from ase.dft.kpoints import resolve_custom_points special_points = dict(special_points) path = resolve_custom_points(path, special_points, self._eps) cell = self.tocell() if transformation is not None: cell = transformation.dot(cell) bandpath = BandPath(cell=cell, path=path, special_points=special_points) return bandpath.interpolate(npoints=npoints, density=density)
@abstractmethod def _cell(self, **kwargs): """Return a Cell object from this Bravais lattice. Arguments are the dictionary of Bravais parameters.""" pass def _special_points(self, **kwargs): """Return the special point coordinates as an npoints x 3 sequence. Subclasses typically return a nested list. Ordering must be same as kpoint labels. Arguments are the dictionary of Bravais parameters and the variant.""" raise NotImplementedError def _variant_name(self, **kwargs): """Return the name (e.g. 'ORCF3') of variant. Arguments will be the dictionary of Bravais parameters.""" raise NotImplementedError def __format__(self, spec): tokens = [] if not spec: spec = '.6g' template = '{}={:%s}' % spec for name in self.parameters: value = self._parameters[name] tokens.append(template.format(name, value)) return '{}({})'.format(self.name, ', '.join(tokens)) def __str__(self): return self.__format__('') def __repr__(self): return self.__format__('.20g')
[docs] def description(self): """Return complete description of lattice and Brillouin zone.""" points = self.get_special_points() labels = self.special_point_names coordstring = '\n'.join([' {:2s} {:7.4f} {:7.4f} {:7.4f}' .format(label, *points[label]) for label in labels]) string = """\ {repr} {variant} Special point coordinates: {special_points} """.format(repr=str(self), variant=self._variant, special_points=coordstring) return string
[docs] @classmethod def type_description(cls): """Return complete description of this Bravais lattice type.""" desc = """\ Lattice name: {name} Long name: {longname} Parameters: {parameters} """.format(**vars(cls)) chunks = [desc] for name in cls.variant_names: var = cls.variants[name] txt = str(var) lines = [' ' + L for L in txt.splitlines()] lines.append('') chunks.extend(lines) return '\n'.join(chunks)
class Variant: variant_desc = """\ Variant name: {name} Special point names: {special_point_names} Default path: {special_path} """ def __init__(self, name, special_point_names, special_path, special_points=None): self.name = name self.special_point_names = special_point_names self.special_path = special_path if special_points is not None: special_points = dict(special_points) for key, value in special_points.items(): special_points[key] = np.array(value) self.special_points = special_points # XXX Should make special_points available as a single array as well # (easier to transform that way) def __str__(self): return self.variant_desc.format(**vars(self)) bravais_names = [] bravais_lattices = {} bravais_classes = {} def bravaisclass(longname, crystal_family, lattice_system, pearson_symbol, parameters, variants, ndim=3): """Decorator for Bravais lattice classes. This sets a number of class variables and processes the information about different variants into a list of Variant objects.""" def decorate(cls): btype = cls.__name__ cls.name = btype cls.longname = longname cls.crystal_family = crystal_family cls.lattice_system = lattice_system cls.pearson_symbol = pearson_symbol cls.parameters = tuple(parameters) cls.variant_names = [] cls.variants = {} cls.ndim = ndim for [name, special_point_names, special_path, special_points] in variants: cls.variant_names.append(name) cls.variants[name] = Variant(name, special_point_names, special_path, special_points) # Register in global list and dictionary bravais_names.append(btype) bravais_lattices[btype] = cls bravais_classes[pearson_symbol] = cls return cls return decorate class UnconventionalLattice(ValueError): pass class Cubic(BravaisLattice): """Abstract class for cubic lattices.""" def __init__(self, a): BravaisLattice.__init__(self, a=a) @bravaisclass('primitive cubic', 'cubic', 'cubic', 'cP', 'a', [['CUB', 'GXRM', 'GXMGRX,MR', sc_special_points['cubic']]]) class CUB(Cubic): def _cell(self, a): return a * np.eye(3) @bravaisclass('face-centred cubic', 'cubic', 'cubic', 'cF', 'a', [['FCC', 'GKLUWX', 'GXWKGLUWLK,UX', sc_special_points['fcc']]]) class FCC(Cubic): def _cell(self, a): return 0.5 * np.array([[0., a, a], [a, 0, a], [a, a, 0]]) @bravaisclass('body-centred cubic', 'cubic', 'cubic', 'cI', 'a', [['BCC', 'GHPN', 'GHNGPH,PN', sc_special_points['bcc']]]) class BCC(Cubic): def _cell(self, a): return 0.5 * np.array([[-a, a, a], [a, -a, a], [a, a, -a]]) @bravaisclass('primitive tetragonal', 'tetragonal', 'tetragonal', 'tP', 'ac', [['TET', 'GAMRXZ', 'GXMGZRAZ,XR,MA', sc_special_points['tetragonal']]]) class TET(BravaisLattice): def __init__(self, a, c): BravaisLattice.__init__(self, a=a, c=c) def _cell(self, a, c): return np.diag(np.array([a, a, c])) # XXX in BCT2 we use S for Sigma. # Also in other places I think @bravaisclass('body-centred tetragonal', 'tetragonal', 'tetragonal', 'tI', 'ac', [['BCT1', 'GMNPXZZ1', 'GXMGZPNZ1M,XP', None], ['BCT2', 'GNPSS1XYY1Z', 'GXYSGZS1NPY1Z,XP', None]]) class BCT(BravaisLattice): def __init__(self, a, c): BravaisLattice.__init__(self, a=a, c=c) def _cell(self, a, c): return 0.5 * np.array([[-a, a, c], [a, -a, c], [a, a, -c]]) def _variant_name(self, a, c): return 'BCT1' if c < a else 'BCT2' def _special_points(self, a, c, variant): a2 = a * a c2 = c * c assert variant.name in self.variants if variant.name == 'BCT1': eta = .25 * (1 + c2 / a2) points = [[0,0,0], [-.5, .5, .5], [0.,.5,0.], [.25, .25, .25], [0.,0.,.5], [eta,eta,-eta], [-eta,1-eta,eta]] else: eta = .25 * (1 + a2 / c2) # Not same eta as BCT1! zeta = 0.5 * a2 / c2 points = [[0.,.0,0.], [0.,.5,0.], [.25,.25,.25], [-eta,eta,eta], [eta,1-eta,-eta], [0.,0.,.5], [-zeta,zeta,.5], [.5,.5,-zeta], [.5,.5,-.5]] return points def check_orc(a, b, c): if not a < b < c: raise UnconventionalLattice('Expected a < b < c, got {}, {}, {}' .format(a, b, c)) class Orthorhombic(BravaisLattice): """Abstract class for orthorhombic types.""" def __init__(self, a, b, c): check_orc(a, b, c) BravaisLattice.__init__(self, a=a, b=b, c=c) @bravaisclass('primitive orthorhombic', 'orthorhombic', 'orthorhombic', 'oP', 'abc', [['ORC', 'GRSTUXYZ', 'GXSYGZURTZ,YT,UX,SR', sc_special_points['orthorhombic']]]) class ORC(Orthorhombic): def _cell(self, a, b, c): return np.diag([a, b, c]).astype(float) @bravaisclass('face-centred orthorhombic', 'orthorhombic', 'orthorhombic', 'oF', 'abc', [['ORCF1', 'GAA1LTXX1YZ', 'GYTZGXA1Y,TX1,XAZ,LG', None], ['ORCF2', 'GCC1DD1LHH1XYZ', 'GYCDXGZD1HC,C1Z,XH1,HY,LG', None], ['ORCF3', 'GAA1LTXX1YZ', 'GYTZGXA1Y,XAZ,LG', None]]) class ORCF(Orthorhombic): def _cell(self, a, b, c): return 0.5 * np.array([[0, b, c], [a, 0, c], [a, b, 0]]) def _special_points(self, a, b, c, variant): a2 = a * a b2 = b * b c2 = c * c xminus = 0.25 * (1 + a2 / b2 - a2 / c2) xplus = 0.25 * (1 + a2 / b2 + a2 / c2) if variant.name == 'ORCF1' or variant.name == 'ORCF3': zeta = xminus eta = xplus points = [[0, 0, 0], [.5, .5 + zeta, zeta], [.5, .5 - zeta, 1 - zeta], [.5, .5, .5], [1., .5, .5], [0., eta, eta], [1., 1 - eta, 1 - eta], [.5, 0, .5], [.5, .5, 0]] else: assert variant.name == 'ORCF2' phi = 0.25 * (1 + c2 / b2 - c2 / a2) delta = 0.25 * (1 + b2 / a2 - b2 / c2) eta = xminus points = [[0,0,0], [.5, .5-eta, 1-eta], [.5, .5+eta, eta], [.5-delta, .5, 1-delta], [.5+delta, .5, delta], [.5, .5, .5], [1-phi, .5-phi, .5], [phi, .5+phi, .5], [0., .5, .5], [.5, 0., .5], [.5, .5, 0.]] return points def _variant_name(self, a, b, c): diff = 1.0 / (a * a) - 1.0 / (b * b) - 1.0 / (c * c) if abs(diff) < self._eps: return 'ORCF3' return 'ORCF1' if diff > 0 else 'ORCF2' @bravaisclass('body-centred orthorhombic', 'orthorhombic', 'orthorhombic', 'oI', 'abc', [['ORCI', 'GLL1L2RSTWXX1YY1Z', 'GXLTWRX1ZGYSW,L1Y,Y1Z', None]]) class ORCI(Orthorhombic): def _cell(self, a, b, c): return 0.5 * np.array([[-a, b, c], [a, -b, c], [a, b, -c]]) def _special_points(self, a, b, c, variant): a2 = a**2 b2 = b**2 c2 = c**2 zeta = .25 * (1 + a2 / c2) eta = .25 * (1 + b2 / c2) delta = .25 * (b2 - a2) / c2 mu = .25 * (a2 + b2) / c2 points = [[0.,0.,0.], [-mu,mu,.5-delta], [mu, -mu, .5+delta], [.5-delta, .5+delta, -mu], [0,.5,0], [.5,0,0], [0.,0.,.5], [.25,.25,.25], [-zeta, zeta, zeta], [zeta, 1 - zeta, -zeta], [eta, -eta, eta], [1 - eta, eta, -eta], [.5,.5,-.5]] return points @bravaisclass('base-centred orthorhombic', 'orthorhombic', 'orthorhombic', 'oC', 'abc', [['ORCC', 'GAA1RSTXX1YZ', 'GXSRAZGYX1A1TY,ZT', None]]) class ORCC(BravaisLattice): def __init__(self, a, b, c): # ORCC is the only ORCx lattice with a<b and not a<b<c if a >= b: raise UnconventionalLattice('Expected a < b, got {}, {}' .format(a, b, c)) BravaisLattice.__init__(self, a=a, b=b, c=c) def _cell(self, a, b, c): return np.array([[0.5 * a, -0.5 * b, 0], [0.5 * a, 0.5 * b, 0], [0, 0, c]]) def _special_points(self, a, b, c, variant): zeta = .25 * (1 + a * a / (b * b)) points = [[0,0,0], [zeta,zeta,.5], [-zeta,1-zeta,.5], [0,.5,.5], [0,.5,0], [-.5,.5,.5], [zeta,zeta,0], [-zeta,1-zeta,0], [-.5,.5,0], [0,0,.5]] return points @bravaisclass('primitive hexagonal', 'hexagonal', 'hexagonal', 'hP', 'ac', [['HEX', 'GMKALH', 'GMKGALHA,LM,KH', sc_special_points['hexagonal']]]) class HEX(BravaisLattice): def __init__(self, a, c): BravaisLattice.__init__(self, a=a, c=c) def _cell(self, a, c): x = 0.5 * np.sqrt(3) return np.array([[0.5 * a, -x * a, 0], [0.5 * a, x * a, 0], [0., 0., c]]) @bravaisclass('primitive rhombohedral', 'hexagonal', 'rhombohedral', 'hR', ('a', 'alpha'), [['RHL1', 'GBB1FLL1PP1P2QXZ', 'GLB1,BZGX,QFP1Z,LP', None], ['RHL2', 'GFLPP1QQ1Z', 'GPZQGFP1Q1LZ', None]]) class RHL(BravaisLattice): def __init__(self, a, alpha): if alpha >= 120: raise UnconventionalLattice('Need alpha < 120 degrees, got {}' .format(alpha)) BravaisLattice.__init__(self, a=a, alpha=alpha) def _cell(self, a, alpha): alpha *= np.pi / 180 acosa = a * np.cos(alpha) acosa2 = a * np.cos(0.5 * alpha) asina2 = a * np.sin(0.5 * alpha) acosfrac = acosa / acosa2 xx = (1 - acosfrac**2) assert xx > 0.0 return np.array([[acosa2, -asina2, 0], [acosa2, asina2, 0], [a * acosfrac, 0, a * xx**0.5]]) def _variant_name(self, a, alpha): return 'RHL1' if alpha < 90 else 'RHL2' def _special_points(self, a, alpha, variant): if variant.name == 'RHL1': cosa = np.cos(alpha * _degrees) eta = (1 + 4 * cosa) / (2 + 4 * cosa) nu = .75 - 0.5 * eta points = [[0,0,0], [eta,.5,1-eta], [.5, 1 - eta, eta - 1], [.5,.5,0], [.5,0,0], [0,0,-.5], [eta,nu,nu], [1-nu,1-nu,1-eta], [nu,nu,eta-1], [1-nu,nu,0], [nu,0,-nu], [.5,.5,.5]] else: eta = 1 / (2 * np.tan(alpha * _degrees / 2)**2) nu = .75 - 0.5 * eta points = [[0,0,0], [.5,-.5,0], [.5,0,0], [1-nu,-nu,1-nu], [nu,nu-1,nu-1], [eta,eta,eta], [1-eta,-eta,-eta], [.5,-.5,.5]] return points def check_mcl(a, b, c, alpha): if not (b <= c and alpha < 90): raise UnconventionalLattice('Expected b <= c, alpha < 90; ' 'got a={}, b={}, c={}, alpha={}' .format(a, b, c, alpha)) @bravaisclass('primitive monoclinic', 'monoclinic', 'monoclinic', 'mP', ('a', 'b', 'c', 'alpha'), [['MCL', 'GACDD1EHH1H2MM1M2XYY1Z', 'GYHCEM1AXH1,MDZ,YD', None]]) class MCL(BravaisLattice): def __init__(self, a, b, c, alpha): check_mcl(a, b, c, alpha) BravaisLattice.__init__(self, a=a, b=b, c=c, alpha=alpha) def _cell(self, a, b, c, alpha): alpha *= _degrees return np.array([[a, 0, 0], [0, b, 0], [0, c * np.cos(alpha), c * np.sin(alpha)]]) def _special_points(self, a, b, c, alpha, variant): cosa = np.cos(alpha * _degrees) eta = (1 - b * cosa / c) / (2 * np.sin(alpha * _degrees)**2) nu = .5 - eta * c * cosa / b points = [[0,0,0], [.5,.5,0], [0,.5,.5], [.5,0,.5], [.5,0,-.5], [.5,.5,.5], [0,eta,1-nu], [0,1-eta,nu], [0,eta,-nu], [.5,eta,1-nu], [.5,1-eta,nu], [.5,eta,-nu], [0,.5,0], [0,0,.5], [0,0,-.5], [.5,0,0]] return points def _variant_name(self, a, b, c, alpha): check_mcl(a, b, c, alpha) return 'MCL' @bravaisclass('base-centred monoclinic', 'monoclinic', 'monoclinic', 'mC', ('a', 'b', 'c', 'alpha'), [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z', 'GYFLI,I1ZF1,YX1,XGN,MG', None], ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z', 'GYFLI,I1ZF1,NGM', None], ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z', 'GYFHZIF1,H1Y1XGN,MG', None], ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z', 'GYFHZI,H1Y1XGN,MG', None], ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z', 'GYFLI,I1ZHF1,H1Y1XGN,MG', None]]) class MCLC(BravaisLattice): def __init__(self, a, b, c, alpha): check_mcl(a, b, c, alpha) BravaisLattice.__init__(self, a=a, b=b, c=c, alpha=alpha) def _cell(self, a, b, c, alpha): alpha *= np.pi / 180 return np.array([[0.5 * a, 0.5 * b, 0], [-0.5 * a, 0.5 * b, 0], [0, c * np.cos(alpha), c * np.sin(alpha)]]) def _variant_name(self, a, b, c, alpha): #from ase.geometry.cell import mclc # okay, this is a bit hacky # We need the same parameters here as when determining the points. # Right now we just repeat the code: check_mcl(a, b, c, alpha) a2 = a * a b2 = b * b cosa = np.cos(alpha * _degrees) sina = np.sin(alpha * _degrees) sina2 = sina**2 cell = self.tocell() lengths_angles = Cell(cell.reciprocal()).cellpar() kgamma = lengths_angles[-1] eps = self._eps # We should not compare angles in degrees versus lengths with # the same precision. if abs(kgamma - 90) < eps: variant = 2 elif kgamma > 90: variant = 1 elif kgamma < 90: num = b * cosa / c + b2 * sina2 / a2 if abs(num - 1) < eps: variant = 4 elif num < 1: variant = 3 else: variant = 5 variant = 'MCLC' + str(variant) return variant def _special_points(self, a, b, c, alpha, variant): variant = int(variant.name[-1]) a2 = a * a b2 = b * b # c2 = c * c cosa = np.cos(alpha * _degrees) sina = np.sin(alpha * _degrees) sina2 = sina**2 if variant == 1 or variant == 2: zeta = (2 - b * cosa / c) / (4 * sina2) eta = 0.5 + 2 * zeta * c * cosa / b psi = .75 - a2 / (4 * b2 * sina * sina) phi = psi + (.75 - psi) * b * cosa / c points = [[0,0,0], [.5,0,0], [0,-.5,0], [1-zeta,1-zeta,1-eta], [zeta,zeta,eta], [-zeta,-zeta,1-eta], [1-zeta,-zeta,1-eta], [phi,1-phi,.5], [1-phi,phi-1,.5], [.5,.5,.5], [.5,0,.5], [1-psi,psi-1,0], [psi,1-psi,0], [psi-1,-psi,0], [.5,.5,0], [-.5,-.5,0], [0,0,.5]] elif variant == 3 or variant == 4: mu = .25 * (1 + b2 / a2) delta = b * c * cosa / (2 * a2) zeta = mu - 0.25 + (1 - b * cosa / c) / (4 * sina2) eta = 0.5 + 2 * zeta * c * cosa / b phi = 1 + zeta - 2 * mu psi = eta - 2 * delta points = [[0,0,0], [1-phi,1-phi,1-psi], [phi,phi-1,psi], [1-phi,-phi,1-psi], [zeta,zeta,eta], [1-zeta,-zeta,1-eta], [-zeta,-zeta,1-eta], [.5,-.5,.5], [.5,0,.5], [.5,0,0], [0,-.5,0], [.5,-.5,0], [mu,mu,delta], [1-mu,-mu,-delta], [-mu,-mu,-delta], [mu,mu-1,delta], [0,0,.5]] elif variant == 5: zeta = .25 * (b2 / a2 + (1 - b * cosa / c) / sina2) eta = 0.5 + 2 * zeta * c * cosa / b mu = .5 * eta + b2 / (4 * a2) - b * c * cosa / (2 * a2) nu = 2 * mu - zeta omega = (4 * nu - 1 - b2 * sina2 / a2) * c / (2 * b * cosa) delta = zeta * c * cosa / b + omega / 2 - .25 rho = 1 - zeta * a2 / b2 points = [[0,0,0], [nu,nu,omega], [1-nu,1-nu,1-omega], [nu,nu-1,omega], [zeta,zeta,eta], [1-zeta,-zeta,1-eta], [-zeta,-zeta,1-eta], [rho,1-rho,.5], [1-rho,rho-1,.5], [.5,.5,.5], [.5,0,.5], [.5,0,0], [0,-.5,0], [.5,-.5,0], [mu,mu,delta], [1-mu,-mu,-delta], [-mu,-mu,-delta], [mu,mu-1,delta], [0,0,.5]] return points tri_angles_explanation = """\ Angles kalpha, kbeta and kgamma of TRI lattice must be 1) all greater \ than 90 degrees with kgamma being the smallest, or 2) all smaller than \ 90 with kgamma being the largest, or 3) kgamma=90 being the \ smallest of the three, or 4) kgamma=90 being the largest of the three. \ Angles of reciprocal lattice are kalpha={}, kbeta={}, kgamma={}. \ If you don't care, please use Cell.fromcellpar() instead.""" # XXX labels, paths, are all the same. @bravaisclass('primitive triclinic', 'triclinic', 'triclinic', 'aP', ('a', 'b', 'c', 'alpha', 'beta', 'gamma'), [['TRI1a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None], ['TRI2a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None], ['TRI1b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None], ['TRI2b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None]]) class TRI(BravaisLattice): def __init__(self, a, b, c, alpha, beta, gamma): BravaisLattice.__init__(self, a=a, b=b, c=c, alpha=alpha, beta=beta, gamma=gamma) def _cell(self, a, b, c, alpha, beta, gamma): alpha, beta, gamma = np.array([alpha, beta, gamma]) singamma = np.sin(gamma * _degrees) cosgamma = np.cos(gamma * _degrees) cosbeta = np.cos(beta * _degrees) cosalpha = np.cos(alpha * _degrees) a3x = c * cosbeta a3y = c / singamma * (cosalpha - cosbeta * cosgamma) a3z = c / singamma * np.sqrt(singamma**2 - cosalpha**2 - cosbeta**2 + 2 * cosalpha * cosbeta * cosgamma) return np.array([[a, 0, 0], [b * cosgamma, b * singamma, 0], [a3x, a3y, a3z]]) def _variant_name(self, a, b, c, alpha, beta, gamma): cell = Cell.new([a, b, c, alpha, beta, gamma]) icellpar = Cell(cell.reciprocal()).cellpar() kangles = kalpha, kbeta, kgamma = icellpar[3:] def raise_unconventional(): raise UnconventionalLattice(tri_angles_explanation .format(*kangles)) eps = self._eps if abs(kgamma - 90) < eps: if kalpha > 90 and kbeta > 90: var = '2a' elif kalpha < 90 and kbeta < 90: var = '2b' else: # Is this possible? Maybe due to epsilon raise_unconventional() elif all(kangles > 90): if kgamma > min(kangles): raise_unconventional() var = '1a' elif all(kangles < 90):# and kgamma > max(kalpha, kbeta): if kgamma < max(kangles): raise_unconventional() var = '1b' else: raise_unconventional() return 'TRI' + var def _special_points(self, a, b, c, alpha, beta, gamma, variant): # (None of the points actually depend on any parameters) # (We should store the points openly on the variant objects) if variant.name == 'TRI1a' or variant.name == 'TRI2a': points = [[0.,0.,0.], [.5,.5,0], [0,.5,.5], [.5,0,.5], [.5,.5,.5], [.5,0,0], [0,.5,0], [0,0,.5]] else: points = [[0,0,0], [.5,-.5,0], [0,0,.5], [-.5,-.5,.5], [0,-.5,.5], [0,-0.5,0], [.5,0,0], [-.5,0,.5]] return points def get_subset_points(names, points): newpoints = {} for name in names: newpoints[name] = points[name] return newpoints @bravaisclass('primitive oblique', 'monoclinic', None, 'mp', ('a', 'b', 'alpha'), [['OBL', 'GYHCH1X', 'GYHCH1XG', None]], ndim=2) class OBL(BravaisLattice): def __init__(self, a, b, alpha, **kwargs): BravaisLattice.__init__(self, a=a, b=b, alpha=alpha, **kwargs) def _cell(self, a, b, alpha): cosa = np.cos(alpha * _degrees) sina = np.sin(alpha * _degrees) return np.array([[a, 0, 0], [b * cosa, b * sina, 0], [0., 0., 0.]]) def _special_points(self, a, b, alpha, variant): # XXX Check me if alpha > 90: _alpha = 180 - alpha a, b = b, a else: _alpha = alpha cosa = np.cos(_alpha * _degrees) eta = (1 - a * cosa / b) / (2 * np.sin(_alpha * _degrees)**2) nu = .5 - eta * b * cosa / a points = [[0, 0, 0], [0, 0.5, 0], [eta, 1 - nu, 0], [.5, .5, 0], [1 - eta, nu, 0], [.5, 0, 0]] if alpha > 90: # Then we map the special points back op = np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]]) points = np.dot(points, op.T) return points @bravaisclass('primitive hexagonal', 'hexagonal', None, 'hp', 'a', [['HEX2D', 'GMK', 'GMKG', get_subset_points('GMK', sc_special_points['hexagonal'])]], ndim=2) class HEX2D(BravaisLattice): def __init__(self, a, **kwargs): BravaisLattice.__init__(self, a=a, **kwargs) def _cell(self, a): x = 0.5 * np.sqrt(3) return np.array([[a, 0, 0], [-0.5 * a, x * a, 0], [0., 0., 0.]]) @bravaisclass('primitive rectangular', 'orthorhombic', None, 'op', 'ab', [['RECT', 'GXSY', 'GXSYGS', get_subset_points('GXSY', sc_special_points['orthorhombic'])]], ndim=2) class RECT(BravaisLattice): def __init__(self, a, b, **kwargs): BravaisLattice.__init__(self, a=a, b=b, **kwargs) def _cell(self, a, b): return np.array([[a, 0, 0], [0, b, 0], [0, 0, 0.]]) @bravaisclass('centred rectangular', 'orthorhombic', None, 'oc', ('a', 'alpha'), [['CRECT', 'GXA1Y', 'GXA1YG', None]], ndim=2) class CRECT(BravaisLattice): def __init__(self, a, alpha, **kwargs): BravaisLattice.__init__(self, a=a, alpha=alpha, **kwargs) def _cell(self, a, alpha): x = np.cos(alpha * _degrees) y = np.sin(alpha * _degrees) return np.array([[a, 0, 0], [a * x, a * y, 0], [0, 0, 0.]]) def _special_points(self, a, alpha, variant): if alpha > 90: _alpha = 180 - alpha else: _alpha = alpha sina2 = np.sin(_alpha / 2 * _degrees)**2 sina = np.sin(_alpha * _degrees)**2 eta = sina2 / sina cosa = np.cos(_alpha * _degrees) xi = eta * cosa points = [[0, 0, 0], [eta, - eta, 0], [0.5 + xi, 0.5 - xi, 0], [0.5, 0.5, 0]] if alpha > 90: # Then we map the special points back op = np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]]) points = np.dot(points, op.T) return points @bravaisclass('primitive square', 'tetragonal', None, 'tp', ('a',), [['SQR', 'GMX', 'MGXM', get_subset_points('GMX', sc_special_points['tetragonal'])]], ndim=2) class SQR(BravaisLattice): def __init__(self, a, **kwargs): BravaisLattice.__init__(self, a=a, **kwargs) def _cell(self, a): return np.array([[a, 0, 0], [0, a, 0], [0, 0, 0.]]) @bravaisclass('primitive line', 'line', None, '?', ('a',), [['LINE', 'GX', 'GX', {'G': [0, 0, 0], 'X': [0.5, 0, 0]}]], ndim=1) class LINE(BravaisLattice): def __init__(self, a, **kwargs): BravaisLattice.__init__(self, a=a, **kwargs) def _cell(self, a): return np.array([[a, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]) def celldiff(cell1, cell2): """Return a unitless measure of the difference between two cells.""" cell1 = Cell.ascell(cell1).complete() cell2 = Cell.ascell(cell2).complete() v1v2 = cell1.volume * cell2.volume if v1v2 == 0: raise ZeroDivisionError('Cell volumes are zero') scale = v1v2**(-1. / 3.) # --> 1/Ang^2 x1 = cell1 @ cell1.T x2 = cell2 @ cell2.T dev = scale * np.abs(x2 - x1).max() return dev def get_lattice_from_canonical_cell(cell, eps=2e-4): """Return a Bravais lattice representing the given cell. This works only for cells that are derived from the standard form (as generated by lat.tocell()) or rotations thereof. If the given cell does not resemble the known form of a Bravais lattice, raise RuntimeError.""" return LatticeChecker(cell, eps).match() def identify_lattice(cell, eps=2e-4, *, pbc=True): """Find Bravais lattice representing this cell. Returns Bravais lattice object representing the cell along with an operation that, applied to the cell, yields the same lengths and angles as the Bravais lattice object.""" pbc = cell.any(1) & pbc2pbc(pbc) npbc = sum(pbc) if npbc == 1: i = np.argmax(pbc) # index of periodic axis a = cell[i, i] if a < 0 or cell[i, [i - 1, i - 2]].any(): raise ValueError('Not a 1-d cell ASE can handle: {cell}.' .format(cell=cell)) if i == 0: op = np.eye(3) elif i == 1: op = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 1]]) else: op = np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]]) return LINE(a), op if npbc == 2: lat, op = get_2d_bravais_lattice(cell, eps, pbc=pbc) return lat, op if npbc != 3: raise ValueError('System must be periodic either ' 'along all three axes, ' 'along two first axes or, ' 'along the thrid axis. ' 'Got pbc={}'.format(pbc)) from ase.geometry.bravais_type_engine import niggli_op_table if cell.rank < 3: raise ValueError('Expected 3 linearly independent cell vectors') rcell, reduction_op = cell.niggli_reduce(eps=eps) # We tabulate the cell's Niggli-mapped versions so we don't need to # redo any work when the same Niggli-operation appears multiple times # in the table: memory = {} # We loop through the most symmetric kinds (CUB etc.) and return # the first one we find: for latname in LatticeChecker.check_order: # There may be multiple Niggli operations that produce valid # lattices, at least for MCL. In that case we will pick the # one whose angle is closest to 90, but it means we cannot # just return the first one we find so we must remember then: matching_lattices = [] for op_key in niggli_op_table[latname]: checker_and_op = memory.get(op_key) if checker_and_op is None: normalization_op = np.array(op_key).reshape(3, 3) candidate = Cell(np.linalg.inv(normalization_op.T) @ rcell) checker = LatticeChecker(candidate, eps=eps) memory[op_key] = (checker, normalization_op) else: checker, normalization_op = checker_and_op lat = checker.query(latname) if lat is not None: op = normalization_op @ np.linalg.inv(reduction_op) matching_lattices.append((lat, op)) # Among any matching lattices, return the one with lowest # orthogonality defect: best = None best_defect = np.inf for lat, op in matching_lattices: cell = lat.tocell() lengths = cell.lengths() defect = np.prod(lengths) / cell.volume if defect < best_defect: best = lat, op best_defect = defect if best is not None: return best class LatticeChecker: # The check order is slightly different than elsewhere listed order # as we need to check HEX/RHL before the ORCx family. check_order = ['CUB', 'FCC', 'BCC', 'TET', 'BCT', 'HEX', 'RHL', 'ORC', 'ORCF', 'ORCI', 'ORCC', 'MCL', 'MCLC', 'TRI'] def __init__(self, cell, eps=2e-4): """Generate Bravais lattices that look (or not) like the given cell. The cell must be reduced to canonical form, i.e., it must be possible to produce a cell with the same lengths and angles by directly through one of the Bravais lattice classes. Generally for internal use (this module). For each of the 14 Bravais lattices, this object can produce a lattice object which represents the same cell, or None if the tolerance eps is not met.""" self.cell = cell self.eps = eps self.cellpar = cell.cellpar() self.lengths = self.A, self.B, self.C = self.cellpar[:3] self.angles = self.cellpar[3:] # Use a 'neutral' length for checking cubic lattices self.A0 = self.lengths.mean() # Vector of the diagonal and then off-diagonal dot products: # [a1 · a1, a2 · a2, a3 · a3, a2 · a3, a3 · a1, a1 · a2] self.prods = (cell @ cell.T).flat[[0, 4, 8, 5, 2, 1]] def _check(self, latcls, *args): if any(arg <= 0 for arg in args): return None try: lat = latcls(*args) except UnconventionalLattice: return None newcell = lat.tocell() err = celldiff(self.cell, newcell) if err < self.eps: return lat def match(self): """Match cell against all lattices, returning most symmetric match. Returns the lattice object. Raises RuntimeError on failure.""" for name in self.check_order: lat = self.query(name) if lat: return lat else: raise RuntimeError('Could not find lattice type for cell ' 'with lengths and angles {}' .format(self.cell.cellpar().tolist())) def query(self, latname): """Match cell against named Bravais lattice. Return lattice object on success, None on failure.""" meth = getattr(self, latname) lat = meth() return lat def CUB(self): # These methods (CUB, FCC, ...) all return a lattice object if # it matches, else None. return self._check(CUB, self.A0) def FCC(self): return self._check(FCC, np.sqrt(2) * self.A0) def BCC(self): return self._check(BCC, 2.0 * self.A0 / np.sqrt(3)) def TET(self): return self._check(TET, self.A, self.C) def _bct_orci_lengths(self): # Coordinate-system independent relation for BCT and ORCI # standard cells: # a1 · a1 + a2 · a3 == a² / 2 # a2 · a2 + a3 · a1 == a² / 2 (BCT) # == b² / 2 (ORCI) # a3 · a3 + a1 · a2 == c² / 2 # We use these to get a, b, and c in those cases. prods = self.prods lengthsqr = 2.0 * (prods[:3] + prods[3:]) if any(lengthsqr < 0): return None return np.sqrt(lengthsqr) def BCT(self): lengths = self._bct_orci_lengths() if lengths is None: return None return self._check(BCT, lengths[0], lengths[2]) def HEX(self): return self._check(HEX, self.A, self.C) def RHL(self): return self._check(RHL, self.A, self.angles[0]) def ORC(self): return self._check(ORC, *self.lengths) def ORCF(self): # ORCF standard cell: # a2 · a3 = a²/4 # a3 · a1 = b²/4 # a1 · a2 = c²/4 prods = self.prods if all(prods[3:] > 0): orcf_abc = 2 * np.sqrt(prods[3:]) return self._check(ORCF, *orcf_abc) def ORCI(self): lengths = self._bct_orci_lengths() if lengths is None: return None return self._check(ORCI, *lengths) def _orcc_ab(self): # ORCC: a1 · a1 + a2 · a3 = a²/2 # a2 · a2 - a2 · a3 = b²/2 prods = self.prods orcc_sqr_ab = np.empty(2) orcc_sqr_ab[0] = 2.0 * (prods[0] + prods[5]) orcc_sqr_ab[1] = 2.0 * (prods[1] - prods[5]) if all(orcc_sqr_ab > 0): return np.sqrt(orcc_sqr_ab) def ORCC(self): orcc_lengths_ab = self._orcc_ab() if orcc_lengths_ab is None: return None return self._check(ORCC, *orcc_lengths_ab, self.C) def MCL(self): return self._check(MCL, *self.lengths, self.angles[0]) def MCLC(self): # MCLC is similar to ORCC: orcc_ab = self._orcc_ab() if orcc_ab is None: return None prods = self.prods C = self.C mclc_a, mclc_b = orcc_ab[::-1] # a, b reversed wrt. ORCC mclc_cosa = 2.0 * prods[3] / (mclc_b * C) if -1 < mclc_cosa < 1: mclc_alpha = np.arccos(mclc_cosa) * 180 / np.pi return self._check(MCLC, mclc_a, mclc_b, C, mclc_alpha) def TRI(self): return self._check(TRI, *self.cellpar) class UnsupportedLattice(ValueError): pass def get_2d_bravais_lattice(origcell, eps=2e-4, *, pbc=True): pbc = origcell.any(1) & pbc2pbc(pbc) if list(pbc) != [1, 1, 0]: raise UnsupportedLattice('Can only get 2D Bravais lattice of cell with ' 'pbc==[1, 1, 0]; but we have {}'.format(pbc)) nonperiodic = pbc.argmin() # Start with op = I ops = [np.eye(3)] for i in range(-1, 1): for j in range(-1, 1): op = [[1, j], [i, 1]] if np.abs(np.linalg.det(op)) > 1e-5: # Only touch periodic dirs: op = np.insert(op, nonperiodic, [0, 0], 0) op = np.insert(op, nonperiodic, ~pbc, 1) ops.append(np.array(op)) def allclose(a, b): return np.allclose(a, b, atol=eps) symrank = 0 for op in ops: cell = Cell(op.dot(origcell)) cellpar = cell.cellpar() angles = cellpar[3:] # Find a, b and gamma gamma = angles[~pbc][0] a, b = cellpar[:3][pbc] anglesm90 = np.abs(angles - 90) # Maximum one angle different from 90 deg in 2d please if np.sum(anglesm90 > eps) > 1: continue all_lengths_equal = abs(a - b) < eps if all_lengths_equal: if allclose(gamma, 90): lat = SQR(a) rank = 5 elif allclose(gamma, 120): lat = HEX2D(a) rank = 4 else: lat = CRECT(a, gamma) rank = 3 else: if allclose(gamma, 90): lat = RECT(a, b) rank = 2 else: lat = OBL(a, b, gamma) rank = 1 op = lat.get_transformation(origcell) if not allclose(np.dot(op, lat.tocell())[pbc][:, pbc], origcell.array[pbc][:, pbc]): msg = ('Cannot recognize cell at all somehow! {}, {}, {}'. format(a, b, gamma)) raise RuntimeError(msg) if rank > symrank: finalop = op symrank = rank finallat = lat return finallat, finalop.T def all_variants(include_blunt_angles=True): """For testing and examples; yield all variants of all lattices.""" a, b, c = 3., 4., 5. alpha = 55.0 yield CUB(a) yield FCC(a) yield BCC(a) yield TET(a, c) bct1 = BCT(2 * a, c) bct2 = BCT(a, c) assert bct1.variant == 'BCT1' assert bct2.variant == 'BCT2' yield bct1 yield bct2 yield ORC(a, b, c) a0 = np.sqrt(1.0 / (1 / b**2 + 1 / c**2)) orcf1 = ORCF(0.5 * a0, b, c) orcf2 = ORCF(1.2 * a0, b, c) orcf3 = ORCF(a0, b, c) assert orcf1.variant == 'ORCF1' assert orcf2.variant == 'ORCF2' assert orcf3.variant == 'ORCF3' yield orcf1 yield orcf2 yield orcf3 yield ORCI(a, b, c) yield ORCC(a, b, c) yield HEX(a, c) rhl1 = RHL(a, alpha=55.0) assert rhl1.variant == 'RHL1' yield rhl1 rhl2 = RHL(a, alpha=105.0) assert rhl2.variant == 'RHL2' yield rhl2 # With these lengths, alpha < 65 (or so) would result in a lattice that # could also be represented with alpha > 65, which is more conventional. yield MCL(a, b, c, alpha=70.0) mclc1 = MCLC(a, b, c, 80) assert mclc1.variant == 'MCLC1' yield mclc1 # mclc2 has same special points as mclc1 mclc3 = MCLC(1.8 * a, b, c * 2, 80) assert mclc3.variant == 'MCLC3' yield mclc3 # mclc4 has same special points as mclc3 # XXX We should add MCLC2 and MCLC4 as well. mclc5 = MCLC(b, b, 1.1 * b, 70) assert mclc5.variant == 'MCLC5' yield mclc5 def get_tri(kcellpar): # We build the TRI lattices from cellpars of reciprocal cell icell = Cell.fromcellpar(kcellpar) cellpar = Cell(4 * icell.reciprocal()).cellpar() return TRI(*cellpar) tri1a = get_tri([1., 1.2, 1.4, 120., 110., 100.]) assert tri1a.variant == 'TRI1a' yield tri1a tri1b = get_tri([1., 1.2, 1.4, 50., 60., 70.]) assert tri1b.variant == 'TRI1b' yield tri1b tri2a = get_tri([1., 1.2, 1.4, 120., 110., 90.]) assert tri2a.variant == 'TRI2a' yield tri2a tri2b = get_tri([1., 1.2, 1.4, 50., 60., 90.]) assert tri2b.variant == 'TRI2b' yield tri2b yield OBL(a, b, alpha=alpha) yield RECT(a, b) yield CRECT(a, alpha=alpha) yield HEX2D(a) yield SQR(a) yield LINE(a) if include_blunt_angles: beta = 110 yield OBL(a, b, alpha=beta) yield CRECT(a, alpha=beta)