from abc import ABC, abstractmethod
from typing import Dict, List
import numpy as np
from ase.cell import Cell
from ase.dft.kpoints import BandPath, parse_path_string, sc_special_points
from ase.utils import pbc2pbc
_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) -> str:
"""Return name of lattice variant.
>>> from ase.lattice import BCT
>>> BCT(3, 5).variant
'BCT2'
"""
return self._variant.name
def __getattr__(self, name: str):
if name in self._parameters:
return self._parameters[name]
return self.__getattribute__(name) # Raises error
[docs] def vars(self) -> Dict[str, float]:
"""Get parameter names and values of this lattice as a dictionary."""
return dict(self._parameters)
[docs] def conventional(self) -> 'BravaisLattice':
"""Get the conventional cell corresponding to this lattice."""
cls = bravais_lattices[self.conventional_cls]
return cls(**self._parameters)
[docs] def tocell(self) -> Cell:
"""Return this lattice as a :class:`~ase.cell.Cell` object."""
cell = self._cell(**self._parameters)
return Cell(cell)
[docs] def cellpar(self) -> np.ndarray:
"""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) -> str:
"""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) -> List[str]:
"""Return all special point names as a list of strings.
>>> from ase.lattice import BCT
>>> 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) -> np.ndarray:
"""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) -> Dict[str, np.ndarray]:
"""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) -> BandPath:
"""Return a :class:`~ase.dft.kpoints.BandPath` for this lattice.
See :meth:`ase.cell.Cell.bandpath` for description of parameters.
>>> from ase.lattice import BCT
>>> 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
path, special_points = resolve_custom_points(path,
special_points,
self._eps)
cell = self.tocell()
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."""
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 = f'{{}}={{:{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) -> str:
return self.__format__('')
def __repr__(self) -> str:
return self.__format__('.20g')
[docs] def description(self) -> str:
"""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) -> str:
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
# Common mappings of primitive to conventional cells:
_identity = np.identity(3, int)
_fcc_map = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
_bcc_map = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]])
class UnconventionalLattice(ValueError):
pass
class Cubic(BravaisLattice):
"""Abstract class for cubic lattices."""
conventional_cls = 'CUB'
def __init__(self, a):
super().__init__(a=a)
@bravaisclass('primitive cubic', 'cubic', 'cubic', 'cP', 'a',
[['CUB', 'GXRM', 'GXMGRX,MR', sc_special_points['cubic']]])
class CUB(Cubic):
conventional_cellmap = _identity
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):
conventional_cellmap = _bcc_map
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):
conventional_cellmap = _fcc_map
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):
conventional_cls = 'TET'
conventional_cellmap = _identity
def __init__(self, a, c):
super().__init__(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):
conventional_cls = 'TET'
conventional_cellmap = _fcc_map
def __init__(self, a, c):
super().__init__(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)
super().__init__(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):
conventional_cls = 'ORC'
conventional_cellmap = _identity
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):
conventional_cls = 'ORC'
conventional_cellmap = _bcc_map
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):
conventional_cls = 'ORC'
conventional_cellmap = _fcc_map
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):
conventional_cls = 'ORC'
conventional_cellmap = np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 1]])
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(f'Expected a < b, got a={a}, b={b}')
super().__init__(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):
conventional_cls = 'HEX'
conventional_cellmap = _identity
def __init__(self, a, c):
super().__init__(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):
conventional_cls = 'RHL'
conventional_cellmap = _identity
def __init__(self, a, alpha):
if alpha >= 120:
raise UnconventionalLattice('Need alpha < 120 degrees, got {}'
.format(alpha))
super().__init__(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 b > c or 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):
conventional_cls = 'MCL'
conventional_cellmap = _identity
def __init__(self, a, b, c, alpha):
check_mcl(a, b, c, alpha)
super().__init__(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):
conventional_cls = 'MCL'
conventional_cellmap = np.array([[1, -1, 0], [1, 1, 0], [0, 0, 1]])
def __init__(self, a, b, c, alpha):
check_mcl(a, b, c, alpha)
super().__init__(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):
conventional_cls = 'TRI'
conventional_cellmap = _identity
def __init__(self, a, b, c, alpha, beta, gamma):
super().__init__(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 = {name: points[name] for name in names}
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):
check_rect(a, b)
if alpha >= 90:
raise UnconventionalLattice(
f'Expected alpha < 90, got alpha={alpha}')
super().__init__(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):
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]]
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):
super().__init__(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.]])
def check_rect(a, b):
if a >= b:
raise UnconventionalLattice(f'Expected a < b, got a={a}, b={b}')
@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):
check_rect(a, b)
super().__init__(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):
# It would probably be better to define the CRECT cell
# by (a, b) rather than (a, alpha). Then we can require a < b
# like in ordinary RECT.
#
# In 3D, all lattices in the same family generally take
# identical parameters.
if alpha >= 90:
raise UnconventionalLattice(
f'Expected alpha < 90. Got alpha={alpha}')
super().__init__(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):
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]]
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):
super().__init__(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):
super().__init__(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 < 1e-10:
# (Proposed cell may be linearly dependent)
return np.inf
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."""
from ase.geometry.bravais_type_engine import niggli_op_table
pbc = cell.any(1) & pbc2pbc(pbc)
npbc = sum(pbc)
cell = cell.uncomplete(pbc)
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_orders[npbc]:
# 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()[pbc]
generalized_volume = cell.complete().volume
defect = np.prod(lengths) / generalized_volume
if defect < best_defect:
best = lat, op
best_defect = defect
if best is not None:
if npbc == 2:
# The 3x3 operation may flip the z axis, but then the x/y
# components are necessarily also left-handed which
# means a defacto left-handed 2D bandpath.
#
# We repair this by applying an operation that unflips the
# z axis and interchanges x/y:
if op[2, 2] < 0:
repair_op = np.array([[0, 1, 0], [1, 0, 0], [0, 0, -1]])
op = repair_op @ op
best = lat, op
return best
raise RuntimeError('Failed to recognize lattice')
class LatticeChecker:
# The check order is slightly different than elsewhere listed order
# as we need to check HEX/RHL before the ORCx family.
check_orders = {
1: ['LINE'],
2: ['SQR', 'RECT', 'HEX2D', 'CRECT', 'OBL'],
3: ['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_orders[self.cell.rank]:
lat = self.query(name)
if lat:
return lat
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 LINE(self):
return self._check(LINE, self.lengths[0])
def SQR(self):
return self._check(SQR, self.lengths[0])
def RECT(self):
return self._check(RECT, *self.lengths[:2])
def CRECT(self):
return self._check(CRECT, self.lengths[0], self.angles[2])
def HEX2D(self):
return self._check(HEX2D, self.lengths[0])
def OBL(self):
return self._check(OBL, *self.lengths[:2], self.angles[2])
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
if mclc_b > C:
# XXX Temporary fix for certain otherwise
# unrecognizable lattices.
#
# This error could happen if the input lattice maps to
# something just outside the domain of conventional
# lattices (less than the tolerance). Our solution is to
# propose a nearby conventional lattice instead, which
# will then be accepted if it's close enough.
mclc_b = 0.5 * (mclc_b + C)
C = mclc_b
return self._check(MCLC, mclc_a, mclc_b, C, mclc_alpha)
def TRI(self):
return self._check(TRI, *self.cellpar)
def all_variants():
"""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
# Choose an OBL lattice that round-trip-converts to itself.
# The default a/b/alpha parameters result in another representation
# of the same lattice.
yield OBL(a=3.0, b=3.35, alpha=77.85)
yield RECT(a, b)
yield CRECT(a, alpha=alpha)
yield HEX2D(a)
yield SQR(a)
yield LINE(a)