Coverage for ase/lattice/__init_

1# fmt: off

3from abc import ABC, abstractmethod

4from dataclasses import dataclass

5from functools import cached_property

6from typing import Dict, List

8import numpy as np

10from ase.cell import Cell

11from ase.dft.kpoints import BandPath, parse_path_string, sc_special_points

12from ase.utils import pbc2pbc

14_degrees = np.pi / 180

17def lattice_attr(name):

18 def getter(self):

19 try:

20 return self._parameters[name]

21 except KeyError:

22 raise AttributeError(name) from None

24 getter.__name__ = name

25 return property(getter)

28class BravaisLattice(ABC):

29 """Represent Bravais lattices and data related to the Brillouin zone.

31 There are 14 3D Bravais classes: CUB, FCC, BCC, ..., and TRI, and

32 five 2D classes.

34 Each class stores basic static information:

36 >>> from ase.lattice import FCC, MCL

37 >>> FCC.name

38 'FCC'

39 >>> FCC.longname

40 'face-centred cubic'

41 >>> FCC.pearson_symbol

42 'cF'

43 >>> MCL.parameters

44 ('a', 'b', 'c', 'alpha')

46 Each class can be instantiated with the specific lattice parameters

47 that apply to that lattice:

49 >>> MCL(3, 4, 5, 80)

50 MCL(a=3, b=4, c=5, alpha=80)

52 """

53 # These parameters can be set by the @bravais decorator for a subclass.

54 # (We could also use metaclasses to do this, but that's more abstract)

55 #

56 # We have to fight mypy so now we actually set initial values that are

57 # not None, unfortunately.

58 name: str = '' # e.g. 'CUB', 'BCT', 'ORCF', ...

59 longname: str = '' # e.g. 'cubic', 'body-centred tetragonal', ...

60 parameters: tuple = tuple() # e.g. ('a', 'c')

61 variants: dict | None = None # e.g. {'BCT1': <variant object>,

62 # 'BCT2': <variant object>}

63 ndim: int = 0

64 conventional_cls: str | None = None

66 def __init__(self, **kwargs):

67 p = {}

68 eps = kwargs.pop('eps', 2e-4)

69 for k, v in kwargs.items():

70 p[k] = float(v)

71 assert set(p) == set(self.parameters)

72 self._parameters = p

73 self._eps = eps

75 if len(self.variants) == 1:

76 # If there's only one it has the same name as the lattice type

77 self._variant = self.variants[self.name]

78 else:

79 name = self._variant_name(**self._parameters)

80 self._variant = self.variants[name]

82 @property

83 def variant(self) -> str:

84 """Return name of lattice variant.

86 >>> from ase.lattice import BCT

88 >>> BCT(3, 5).variant

89 'BCT2'

90 """

91 return self._variant.name

93 def vars(self) -> Dict[str, float]:

94 """Get parameter names and values of this lattice as a dictionary."""

95 return dict(self._parameters)

97 def conventional(self) -> 'BravaisLattice':

98 """Get the conventional cell corresponding to this lattice."""

99 cls = bravais_lattices[self.conventional_cls]

100 return cls(**self._parameters)

101

102 def tocell(self) -> Cell:

103 """Return this lattice as a :class:`~ase.cell.Cell` object."""

104 cell = self._cell(**self._parameters)

105 return Cell(cell)

106

107 def cellpar(self) -> np.ndarray:

108 """Get cell lengths and angles as array of length 6.

109

110 See :func:`ase.geometry.Cell.cellpar`."""

111 # (Just a brute-force implementation)

112 cell = self.tocell()

113 return cell.cellpar()

114

115 @property

116 def special_path(self) -> str:

117 """Get default special k-point path for this lattice as a string.

118

119 >>> BCT(3, 5).special_path

120 'GXYSGZS1NPY1Z,XP'

121 """

122 return self._variant.special_path

123

124 @property

125 def special_point_names(self) -> List[str]:

126 """Return all special point names as a list of strings.

127

128 >>> from ase.lattice import BCT

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) -> 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 >>> from ase.lattice import BCT

181

182 >>> BCT(3, 5).bandpath()

183 BandPath(path='GXYSGZS1NPY1Z,XP', cell=[3x3], \

184special_points={GNPSS1XYY1Z}, kpts=[51x3])

185

186 .. note:: This produces the standard band path following AFlow

187 conventions. If your cell does not follow this convention,

188 you will need :meth:`ase.cell.Cell.bandpath` instead or

189 the kpoints may not correspond to your particular cell.

190

191 """

192 if special_points is None:

193 special_points = self.get_special_points()

194

195 if path is None:

196 path = self._variant.special_path

197 elif not isinstance(path, str):

198 from ase.dft.kpoints import resolve_custom_points

199 path, special_points = resolve_custom_points(path,

200 special_points,

201 self._eps)

202

203 cell = self.tocell()

204 bandpath = BandPath(cell=cell, path=path,

205 special_points=special_points)

206 return bandpath.interpolate(npoints=npoints, density=density)

207

208 @abstractmethod

209 def _cell(self, **kwargs):

210 """Return a Cell object from this Bravais lattice.

211

212 Arguments are the dictionary of Bravais parameters."""

213

214 def _special_points(self, **kwargs):

215 """Return the special point coordinates as an npoints x 3 sequence.

216

217 Subclasses typically return a nested list.

218

219 Ordering must be same as kpoint labels.

220

221 Arguments are the dictionary of Bravais parameters and the variant."""

222 raise NotImplementedError

223

224 def _variant_name(self, **kwargs):

225 """Return the name (e.g. 'ORCF3') of variant.

226

227 Arguments will be the dictionary of Bravais parameters."""

228 raise NotImplementedError

229

230 def __format__(self, spec):

231 tokens = []

232 if not spec:

233 spec = '.6g'

234 template = f'{{}}={{:{spec}}}'

235

236 for name in self.parameters:

237 value = self._parameters[name]

238 tokens.append(template.format(name, value))

239 return '{}({})'.format(self.name, ', '.join(tokens))

240

241 def __str__(self) -> str:

242 return self.__format__('')

243

244 def __repr__(self) -> str:

245 return self.__format__('.20g')

246

247 def description(self) -> str:

248 """Return complete description of lattice and Brillouin zone."""

249 points = self.get_special_points()

250 labels = self.special_point_names

251

252 coordstring = '\n'.join([' {:2s} {:7.4f} {:7.4f} {:7.4f}'

253 .format(label, *points[label])

254 for label in labels])

255

256 string = """\

257{repr}

258 {variant}

259 Special point coordinates:

260{special_points}

261""".format(repr=str(self),

262 variant=self._variant,

263 special_points=coordstring)

264 return string

265

266 @classmethod

267 def type_description(cls):

268 """Return complete description of this Bravais lattice type."""

269 desc = """\

270Lattice name: {name}

271 Long name: {longname}

272 Parameters: {parameters}

273""".format(**vars(cls))

274

275 chunks = [desc]

276 for name in cls.variant_names:

277 var = cls.variants[name]

278 txt = str(var)

279 lines = [' ' + L for L in txt.splitlines()]

280 lines.append('')

281 chunks.extend(lines)

282

283 return '\n'.join(chunks)

284

285

286class Variant:

287 variant_desc = """\

288Variant name: {name}

289 Special point names: {special_point_names}

290 Default path: {special_path}

291"""

292

293 def __init__(self, name, special_point_names, special_path,

294 special_points=None):

295 self.name = name

296 self.special_point_names = special_point_names

297 self.special_path = special_path

298 if special_points is not None:

299 special_points = dict(special_points)

300 for key, value in special_points.items():

301 special_points[key] = np.array(value)

302 self.special_points = special_points

303 # XXX Should make special_points available as a single array as well

304 # (easier to transform that way)

305

306 def __str__(self) -> str:

307 return self.variant_desc.format(**vars(self))

308

309

310bravais_names = []

311bravais_lattices = {}

312bravais_classes = {}

313

314

315def bravaisclass(longname, crystal_family, lattice_system, pearson_symbol,

316 parameters, variants, ndim=3):

317 """Decorator for Bravais lattice classes.

318

319 This sets a number of class variables and processes the information

320 about different variants into a list of Variant objects."""

321

322 def decorate(cls):

323 btype = cls.__name__

324 cls.name = btype

325 cls.longname = longname

326 cls.crystal_family = crystal_family

327 cls.lattice_system = lattice_system

328 cls.pearson_symbol = pearson_symbol

329 cls.parameters = tuple(parameters)

330 cls.variant_names = []

331 cls.variants = {}

332 cls.ndim = ndim

333

334 for parameter in parameters:

335 setattr(cls, parameter, lattice_attr(parameter))

336

337 for [name, special_point_names, special_path,

338 special_points] in variants:

339 cls.variant_names.append(name)

340 cls.variants[name] = Variant(name, special_point_names,

341 special_path, special_points)

342

343 # Register in global list and dictionary

344 bravais_names.append(btype)

345 bravais_lattices[btype] = cls

346 bravais_classes[pearson_symbol] = cls

347 return cls

348

349 return decorate

350

351

352# Common mappings of primitive to conventional cells:

353_identity = np.identity(3, int)

354_fcc_map = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])

355_bcc_map = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]])

356

357

358class UnconventionalLattice(ValueError):

359 pass

360

361

362class Cubic(BravaisLattice):

363 """Abstract class for cubic lattices."""

364 conventional_cls = 'CUB'

365

366 def __init__(self, a):

367 super().__init__(a=a)

368

369

370@bravaisclass('primitive cubic', 'cubic', 'cubic', 'cP', 'a',

371 [['CUB', 'GXRM', 'GXMGRX,MR', sc_special_points['cubic']]])

372class CUB(Cubic):

373 conventional_cellmap = _identity

374

375 def _cell(self, a):

376 return a * np.eye(3)

377

378

379@bravaisclass('face-centred cubic', 'cubic', 'cubic', 'cF', 'a',

380 [['FCC', 'GKLUWX', 'GXWKGLUWLK,UX', sc_special_points['fcc']]])

381class FCC(Cubic):

382 conventional_cellmap = _bcc_map

383

384 def _cell(self, a):

385 return 0.5 * np.array([[0., a, a], [a, 0, a], [a, a, 0]])

386

387

388@bravaisclass('body-centred cubic', 'cubic', 'cubic', 'cI', 'a',

389 [['BCC', 'GHPN', 'GHNGPH,PN', sc_special_points['bcc']]])

390class BCC(Cubic):

391 conventional_cellmap = _fcc_map

392

393 def _cell(self, a):

394 return 0.5 * np.array([[-a, a, a], [a, -a, a], [a, a, -a]])

395

396

397@bravaisclass('primitive tetragonal', 'tetragonal', 'tetragonal', 'tP', 'ac',

398 [['TET', 'GAMRXZ', 'GXMGZRAZ,XR,MA',

399 sc_special_points['tetragonal']]])

400class TET(BravaisLattice):

401 conventional_cls = 'TET'

402 conventional_cellmap = _identity

403

404 def __init__(self, a, c):

405 super().__init__(a=a, c=c)

406

407 def _cell(self, a, c):

408 return np.diag(np.array([a, a, c]))

409

410

411# XXX in BCT2 we use S for Sigma.

412# Also in other places I think

413@bravaisclass('body-centred tetragonal', 'tetragonal', 'tetragonal', 'tI',

414 'ac',

415 [['BCT1', 'GMNPXZZ1', 'GXMGZPNZ1M,XP', None],

416 ['BCT2', 'GNPSS1XYY1Z', 'GXYSGZS1NPY1Z,XP', None]])

417class BCT(BravaisLattice):

418 conventional_cls = 'TET'

419 conventional_cellmap = _fcc_map

420

421 def __init__(self, a, c):

422 super().__init__(a=a, c=c)

423

424 def _cell(self, a, c):

425 return 0.5 * np.array([[-a, a, c], [a, -a, c], [a, a, -c]])

426

427 def _variant_name(self, a, c):

428 return 'BCT1' if c < a else 'BCT2'

429

430 def _special_points(self, a, c, variant):

431 a2 = a * a

432 c2 = c * c

433

434 assert variant.name in self.variants

435

436 if variant.name == 'BCT1':

437 eta = .25 * (1 + c2 / a2)

438 points = [[0, 0, 0],

439 [-.5, .5, .5],

440 [0., .5, 0.],

441 [.25, .25, .25],

442 [0., 0., .5],

443 [eta, eta, -eta],

444 [-eta, 1 - eta, eta]]

445 else:

446 eta = .25 * (1 + a2 / c2) # Not same eta as BCT1!

447 zeta = 0.5 * a2 / c2

448 points = [[0., .0, 0.],

449 [0., .5, 0.],

450 [.25, .25, .25],

451 [-eta, eta, eta],

452 [eta, 1 - eta, -eta],

453 [0., 0., .5],

454 [-zeta, zeta, .5],

455 [.5, .5, -zeta],

456 [.5, .5, -.5]]

457 return points

458

459

460def check_orc(a, b, c):

461 if not a < b < c:

462 raise UnconventionalLattice('Expected a < b < c, got {}, {}, {}'

463 .format(a, b, c))

464

465

466class Orthorhombic(BravaisLattice):

467 """Abstract class for orthorhombic types."""

468

469 def __init__(self, a, b, c):

470 check_orc(a, b, c)

471 super().__init__(a=a, b=b, c=c)

472

473

474@bravaisclass('primitive orthorhombic', 'orthorhombic', 'orthorhombic', 'oP',

475 'abc',

476 [['ORC', 'GRSTUXYZ', 'GXSYGZURTZ,YT,UX,SR',

477 sc_special_points['orthorhombic']]])

478class ORC(Orthorhombic):

479 conventional_cls = 'ORC'

480 conventional_cellmap = _identity

481

482 def _cell(self, a, b, c):

483 return np.diag([a, b, c]).astype(float)

484

485

486@bravaisclass('face-centred orthorhombic', 'orthorhombic', 'orthorhombic',

487 'oF', 'abc',

488 [['ORCF1', 'GAA1LTXX1YZ', 'GYTZGXA1Y,TX1,XAZ,LG', None],

489 ['ORCF2', 'GCC1DD1LHH1XYZ', 'GYCDXGZD1HC,C1Z,XH1,HY,LG', None],

490 ['ORCF3', 'GAA1LTXX1YZ', 'GYTZGXA1Y,XAZ,LG', None]])

491class ORCF(Orthorhombic):

492 conventional_cls = 'ORC'

493 conventional_cellmap = _bcc_map

494

495 def _cell(self, a, b, c):

496 return 0.5 * np.array([[0, b, c], [a, 0, c], [a, b, 0]])

497

498 def _special_points(self, a, b, c, variant):

499 a2 = a * a

500 b2 = b * b

501 c2 = c * c

502 xminus = 0.25 * (1 + a2 / b2 - a2 / c2)

503 xplus = 0.25 * (1 + a2 / b2 + a2 / c2)

504

505 if variant.name == 'ORCF1' or variant.name == 'ORCF3':

506 zeta = xminus

507 eta = xplus

508

509 points = [[0, 0, 0],

510 [.5, .5 + zeta, zeta],

511 [.5, .5 - zeta, 1 - zeta],

512 [.5, .5, .5],

513 [1., .5, .5],

514 [0., eta, eta],

515 [1., 1 - eta, 1 - eta],

516 [.5, 0, .5],

517 [.5, .5, 0]]

518 else:

519 assert variant.name == 'ORCF2'

520 phi = 0.25 * (1 + c2 / b2 - c2 / a2)

521 delta = 0.25 * (1 + b2 / a2 - b2 / c2)

522 eta = xminus

523

524 points = [[0, 0, 0],

525 [.5, .5 - eta, 1 - eta],

526 [.5, .5 + eta, eta],

527 [.5 - delta, .5, 1 - delta],

528 [.5 + delta, .5, delta],

529 [.5, .5, .5],

530 [1 - phi, .5 - phi, .5],

531 [phi, .5 + phi, .5],

532 [0., .5, .5],

533 [.5, 0., .5],

534 [.5, .5, 0.]]

535

536 return points

537

538 def _variant_name(self, a, b, c):

539 diff = 1.0 / (a * a) - 1.0 / (b * b) - 1.0 / (c * c)

540 if abs(diff) < self._eps:

541 return 'ORCF3'

542 return 'ORCF1' if diff > 0 else 'ORCF2'

543

544

545@bravaisclass('body-centred orthorhombic', 'orthorhombic', 'orthorhombic',

546 'oI', 'abc',

547 [['ORCI', 'GLL1L2RSTWXX1YY1Z', 'GXLTWRX1ZGYSW,L1Y,Y1Z', None]])

548class ORCI(Orthorhombic):

549 conventional_cls = 'ORC'

550 conventional_cellmap = _fcc_map

551

552 def _cell(self, a, b, c):

553 return 0.5 * np.array([[-a, b, c], [a, -b, c], [a, b, -c]])

554

555 def _special_points(self, a, b, c, variant):

556 a2 = a**2

557 b2 = b**2

558 c2 = c**2

559

560 zeta = .25 * (1 + a2 / c2)

561 eta = .25 * (1 + b2 / c2)

562 delta = .25 * (b2 - a2) / c2

563 mu = .25 * (a2 + b2) / c2

564

565 points = [[0., 0., 0.],

566 [-mu, mu, .5 - delta],

567 [mu, -mu, .5 + delta],

568 [.5 - delta, .5 + delta, -mu],

569 [0, .5, 0],

570 [.5, 0, 0],

571 [0., 0., .5],

572 [.25, .25, .25],

573 [-zeta, zeta, zeta],

574 [zeta, 1 - zeta, -zeta],

575 [eta, -eta, eta],

576 [1 - eta, eta, -eta],

577 [.5, .5, -.5]]

578 return points

579

580

581@bravaisclass('base-centred orthorhombic', 'orthorhombic', 'orthorhombic',

582 'oC', 'abc',

583 [['ORCC', 'GAA1RSTXX1YZ', 'GXSRAZGYX1A1TY,ZT', None]])

584class ORCC(BravaisLattice):

585 conventional_cls = 'ORC'

586 conventional_cellmap = np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 1]])

587

588 def __init__(self, a, b, c):

589 # ORCC is the only ORCx lattice with a<b and not a<b<c

590 if a >= b:

591 raise UnconventionalLattice(f'Expected a < b, got a={a}, b={b}')

592 super().__init__(a=a, b=b, c=c)

593

594 def _cell(self, a, b, c):

595 return np.array([[0.5 * a, -0.5 * b, 0], [0.5 * a, 0.5 * b, 0],

596 [0, 0, c]])

597

598 def _special_points(self, a, b, c, variant):

599 zeta = .25 * (1 + a * a / (b * b))

600 points = [[0, 0, 0],

601 [zeta, zeta, .5],

602 [-zeta, 1 - zeta, .5],

603 [0, .5, .5],

604 [0, .5, 0],

605 [-.5, .5, .5],

606 [zeta, zeta, 0],

607 [-zeta, 1 - zeta, 0],

608 [-.5, .5, 0],

609 [0, 0, .5]]

610 return points

611

612

613@bravaisclass('primitive hexagonal', 'hexagonal', 'hexagonal', 'hP',

614 'ac',

615 [['HEX', 'GMKALH', 'GMKGALHA,LM,KH',

616 sc_special_points['hexagonal']]])

617class HEX(BravaisLattice):

618 conventional_cls = 'HEX'

619 conventional_cellmap = _identity

620

621 def __init__(self, a, c):

622 super().__init__(a=a, c=c)

623

624 def _cell(self, a, c):

625 x = 0.5 * np.sqrt(3)

626 return np.array([[0.5 * a, -x * a, 0], [0.5 * a, x * a, 0],

627 [0., 0., c]])

628

629

630@bravaisclass('primitive rhombohedral', 'hexagonal', 'rhombohedral', 'hR',

631 ('a', 'alpha'),

632 [['RHL1', 'GBB1FLL1PP1P2QXZ', 'GLB1,BZGX,QFP1Z,LP', None],

633 ['RHL2', 'GFLPP1QQ1Z', 'GPZQGFP1Q1LZ', None]])

634class RHL(BravaisLattice):

635 conventional_cls = 'RHL'

636 conventional_cellmap = _identity

637

638 def __init__(self, a, alpha):

639 if alpha >= 120:

640 raise UnconventionalLattice('Need alpha < 120 degrees, got {}'

641 .format(alpha))

642 super().__init__(a=a, alpha=alpha)

643

644 def _cell(self, a, alpha):

645 alpha *= np.pi / 180

646 acosa = a * np.cos(alpha)

647 acosa2 = a * np.cos(0.5 * alpha)

648 asina2 = a * np.sin(0.5 * alpha)

649 acosfrac = acosa / acosa2

650 xx = (1 - acosfrac**2)

651 assert xx > 0.0

652 return np.array([[acosa2, -asina2, 0], [acosa2, asina2, 0],

653 [a * acosfrac, 0, a * xx**0.5]])

654

655 def _variant_name(self, a, alpha):

656 return 'RHL1' if alpha < 90 else 'RHL2'

657

658 def _special_points(self, a, alpha, variant):

659 if variant.name == 'RHL1':

660 cosa = np.cos(alpha * _degrees)

661 eta = (1 + 4 * cosa) / (2 + 4 * cosa)

662 nu = .75 - 0.5 * eta

663 points = [[0, 0, 0],

664 [eta, .5, 1 - eta],

665 [.5, 1 - eta, eta - 1],

666 [.5, .5, 0],

667 [.5, 0, 0],

668 [0, 0, -.5],

669 [eta, nu, nu],

670 [1 - nu, 1 - nu, 1 - eta],

671 [nu, nu, eta - 1],

672 [1 - nu, nu, 0],

673 [nu, 0, -nu],

674 [.5, .5, .5]]

675 else:

676 eta = 1 / (2 * np.tan(alpha * _degrees / 2)**2)

677 nu = .75 - 0.5 * eta

678 points = [[0, 0, 0],

679 [.5, -.5, 0],

680 [.5, 0, 0],

681 [1 - nu, -nu, 1 - nu],

682 [nu, nu - 1, nu - 1],

683 [eta, eta, eta],

684 [1 - eta, -eta, -eta],

685 [.5, -.5, .5]]

686 return points

687

688

689def check_mcl(a, b, c, alpha):

690 if b > c or alpha >= 90:

691 raise UnconventionalLattice('Expected b <= c, alpha < 90; '

692 'got a={}, b={}, c={}, alpha={}'

693 .format(a, b, c, alpha))

694

695

696@bravaisclass('primitive monoclinic', 'monoclinic', 'monoclinic', 'mP',

697 ('a', 'b', 'c', 'alpha'),

698 [['MCL', 'GACDD1EHH1H2MM1M2XYY1Z', 'GYHCEM1AXH1,MDZ,YD', None]])

699class MCL(BravaisLattice):

700 conventional_cls = 'MCL'

701 conventional_cellmap = _identity

702

703 def __init__(self, a, b, c, alpha):

704 check_mcl(a, b, c, alpha)

705 super().__init__(a=a, b=b, c=c, alpha=alpha)

706

707 def _cell(self, a, b, c, alpha):

708 alpha *= _degrees

709 return np.array([[a, 0, 0], [0, b, 0],

710 [0, c * np.cos(alpha), c * np.sin(alpha)]])

711

712 def _special_points(self, a, b, c, alpha, variant):

713 cosa = np.cos(alpha * _degrees)

714 eta = (1 - b * cosa / c) / (2 * np.sin(alpha * _degrees)**2)

715 nu = .5 - eta * c * cosa / b

716

717 points = [[0, 0, 0],

718 [.5, .5, 0],

719 [0, .5, .5],

720 [.5, 0, .5],

721 [.5, 0, -.5],

722 [.5, .5, .5],

723 [0, eta, 1 - nu],

724 [0, 1 - eta, nu],

725 [0, eta, -nu],

726 [.5, eta, 1 - nu],

727 [.5, 1 - eta, nu],

728 [.5, eta, -nu],

729 [0, .5, 0],

730 [0, 0, .5],

731 [0, 0, -.5],

732 [.5, 0, 0]]

733 return points

734

735 def _variant_name(self, a, b, c, alpha):

736 check_mcl(a, b, c, alpha)

737 return 'MCL'

738

739

740@bravaisclass('base-centred monoclinic', 'monoclinic', 'monoclinic', 'mC',

741 ('a', 'b', 'c', 'alpha'),

742 [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

743 'GYFLI,I1ZF1,YX1,XGN,MG', None],

744 ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

745 'GYFLI,I1ZF1,NGM', None],

746 ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

747 'GYFHZIF1,H1Y1XGN,MG', None],

748 ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

749 'GYFHZI,H1Y1XGN,MG', None],

750 ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z',

751 'GYFLI,I1ZHF1,H1Y1XGN,MG', None]])

752class MCLC(BravaisLattice):

753 conventional_cls = 'MCL'

754 conventional_cellmap = np.array([[1, -1, 0], [1, 1, 0], [0, 0, 1]])

755

756 def __init__(self, a, b, c, alpha):

757 check_mcl(a, b, c, alpha)

758 super().__init__(a=a, b=b, c=c, alpha=alpha)

759

760 def _cell(self, a, b, c, alpha):

761 alpha *= np.pi / 180

762 return np.array([[0.5 * a, 0.5 * b, 0], [-0.5 * a, 0.5 * b, 0],

763 [0, c * np.cos(alpha), c * np.sin(alpha)]])

764

765 def _variant_name(self, a, b, c, alpha):

766 # from ase.geometry.cell import mclc

767 # okay, this is a bit hacky

768

769 # We need the same parameters here as when determining the points.

770 # Right now we just repeat the code:

771 check_mcl(a, b, c, alpha)

772

773 a2 = a * a

774 b2 = b * b

775 cosa = np.cos(alpha * _degrees)

776 sina = np.sin(alpha * _degrees)

777 sina2 = sina**2

778

779 cell = self.tocell()

780 lengths_angles = Cell(cell.reciprocal()).cellpar()

781

782 kgamma = lengths_angles[-1]

783

784 eps = self._eps

785 # We should not compare angles in degrees versus lengths with

786 # the same precision.

787 if abs(kgamma - 90) < eps:

788 variant = 2

789 elif kgamma > 90:

790 variant = 1

791 elif kgamma < 90:

792 num = b * cosa / c + b2 * sina2 / a2

793 if abs(num - 1) < eps:

794 variant = 4

795 elif num < 1:

796 variant = 3

797 else:

798 variant = 5

799 variant = 'MCLC' + str(variant)

800 return variant

801

802 def _special_points(self, a, b, c, alpha, variant):

803 variant = int(variant.name[-1])

804

805 a2 = a * a

806 b2 = b * b

807 # c2 = c * c

808 cosa = np.cos(alpha * _degrees)

809 sina = np.sin(alpha * _degrees)

810 sina2 = sina**2

811

812 if variant == 1 or variant == 2:

813 zeta = (2 - b * cosa / c) / (4 * sina2)

814 eta = 0.5 + 2 * zeta * c * cosa / b

815 psi = .75 - a2 / (4 * b2 * sina * sina)

816 phi = psi + (.75 - psi) * b * cosa / c

817

818 points = [[0, 0, 0],

819 [.5, 0, 0],

820 [0, -.5, 0],

821 [1 - zeta, 1 - zeta, 1 - eta],

822 [zeta, zeta, eta],

823 [-zeta, -zeta, 1 - eta],

824 [1 - zeta, -zeta, 1 - eta],

825 [phi, 1 - phi, .5],

826 [1 - phi, phi - 1, .5],

827 [.5, .5, .5],

828 [.5, 0, .5],

829 [1 - psi, psi - 1, 0],

830 [psi, 1 - psi, 0],

831 [psi - 1, -psi, 0],

832 [.5, .5, 0],

833 [-.5, -.5, 0],

834 [0, 0, .5]]

835 elif variant == 3 or variant == 4:

836 mu = .25 * (1 + b2 / a2)

837 delta = b * c * cosa / (2 * a2)

838 zeta = mu - 0.25 + (1 - b * cosa / c) / (4 * sina2)

839 eta = 0.5 + 2 * zeta * c * cosa / b

840 phi = 1 + zeta - 2 * mu

841 psi = eta - 2 * delta

842

843 points = [[0, 0, 0],

844 [1 - phi, 1 - phi, 1 - psi],

845 [phi, phi - 1, psi],

846 [1 - phi, -phi, 1 - psi],

847 [zeta, zeta, eta],

848 [1 - zeta, -zeta, 1 - eta],

849 [-zeta, -zeta, 1 - eta],

850 [.5, -.5, .5],

851 [.5, 0, .5],

852 [.5, 0, 0],

853 [0, -.5, 0],

854 [.5, -.5, 0],

855 [mu, mu, delta],

856 [1 - mu, -mu, -delta],

857 [-mu, -mu, -delta],

858 [mu, mu - 1, delta],

859 [0, 0, .5]]

860 elif variant == 5:

861 zeta = .25 * (b2 / a2 + (1 - b * cosa / c) / sina2)

862 eta = 0.5 + 2 * zeta * c * cosa / b

863 mu = .5 * eta + b2 / (4 * a2) - b * c * cosa / (2 * a2)

864 nu = 2 * mu - zeta

865 omega = (4 * nu - 1 - b2 * sina2 / a2) * c / (2 * b * cosa)

866 delta = zeta * c * cosa / b + omega / 2 - .25

867 rho = 1 - zeta * a2 / b2

868

869 points = [[0, 0, 0],

870 [nu, nu, omega],

871 [1 - nu, 1 - nu, 1 - omega],

872 [nu, nu - 1, omega],

873 [zeta, zeta, eta],

874 [1 - zeta, -zeta, 1 - eta],

875 [-zeta, -zeta, 1 - eta],

876 [rho, 1 - rho, .5],

877 [1 - rho, rho - 1, .5],

878 [.5, .5, .5],

879 [.5, 0, .5],

880 [.5, 0, 0],

881 [0, -.5, 0],

882 [.5, -.5, 0],

883 [mu, mu, delta],

884 [1 - mu, -mu, -delta],

885 [-mu, -mu, -delta],

886 [mu, mu - 1, delta],

887 [0, 0, .5]]

888

889 return points

890

891

892tri_angles_explanation = """\

893Angles kalpha, kbeta and kgamma of TRI lattice must be 1) all greater \

894than 90 degrees with kgamma being the smallest, or 2) all smaller than \

89590 with kgamma being the largest, or 3) kgamma=90 being the \

896smallest of the three, or 4) kgamma=90 being the largest of the three. \

897Angles of reciprocal lattice are kalpha={}, kbeta={}, kgamma={}. \

898If you don't care, please use Cell.fromcellpar() instead."""

899

900# XXX labels, paths, are all the same.

901

902

903@bravaisclass('primitive triclinic', 'triclinic', 'triclinic', 'aP',

904 ('a', 'b', 'c', 'alpha', 'beta', 'gamma'),

905 [['TRI1a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],

906 ['TRI2a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],

907 ['TRI1b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],

908 ['TRI2b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None]])

909class TRI(BravaisLattice):

910 conventional_cls = 'TRI'

911 conventional_cellmap = _identity

912

913 def __init__(self, a, b, c, alpha, beta, gamma):

914 super().__init__(a=a, b=b, c=c, alpha=alpha, beta=beta,

915 gamma=gamma)

916

917 def _cell(self, a, b, c, alpha, beta, gamma):

918 alpha, beta, gamma = np.array([alpha, beta, gamma])

919 singamma = np.sin(gamma * _degrees)

920 cosgamma = np.cos(gamma * _degrees)

921 cosbeta = np.cos(beta * _degrees)

922 cosalpha = np.cos(alpha * _degrees)

923 a3x = c * cosbeta

924 a3y = c / singamma * (cosalpha - cosbeta * cosgamma)

925 a3z = c / singamma * np.sqrt(singamma**2 - cosalpha**2 - cosbeta**2

926 + 2 * cosalpha * cosbeta * cosgamma)

927 return np.array([[a, 0, 0], [b * cosgamma, b * singamma, 0],

928 [a3x, a3y, a3z]])

929

930 def _variant_name(self, a, b, c, alpha, beta, gamma):

931 cell = Cell.new([a, b, c, alpha, beta, gamma])

932 icellpar = Cell(cell.reciprocal()).cellpar()

933 kangles = kalpha, kbeta, kgamma = icellpar[3:]

934

935 def raise_unconventional():

936 raise UnconventionalLattice(tri_angles_explanation

937 .format(*kangles))

938

939 eps = self._eps

940 if abs(kgamma - 90) < eps:

941 if kalpha > 90 and kbeta > 90:

942 var = '2a'

943 elif kalpha < 90 and kbeta < 90:

944 var = '2b'

945 else:

946 # Is this possible? Maybe due to epsilon

947 raise_unconventional()

948 elif all(kangles > 90):

949 if kgamma > min(kangles):

950 raise_unconventional()

951 var = '1a'

952 elif all(kangles < 90): # and kgamma > max(kalpha, kbeta):

953 if kgamma < max(kangles):

954 raise_unconventional()

955 var = '1b'

956 else:

957 raise_unconventional()

958

959 return 'TRI' + var

960

961 def _special_points(self, a, b, c, alpha, beta, gamma, variant):

962 # (None of the points actually depend on any parameters)

963 # (We should store the points openly on the variant objects)

964 if variant.name == 'TRI1a' or variant.name == 'TRI2a':

965 points = [[0., 0., 0.],

966 [.5, .5, 0],

967 [0, .5, .5],

968 [.5, 0, .5],

969 [.5, .5, .5],

970 [.5, 0, 0],

971 [0, .5, 0],

972 [0, 0, .5]]

973 else:

974 points = [[0, 0, 0],

975 [.5, -.5, 0],

976 [0, 0, .5],

977 [-.5, -.5, .5],

978 [0, -.5, .5],

979 [0, -0.5, 0],

980 [.5, 0, 0],

981 [-.5, 0, .5]]

982

983 return points

984

985

986def get_subset_points(names, points):

987 newpoints = {name: points[name] for name in names}

988 return newpoints

989

990

991@bravaisclass('primitive oblique', 'monoclinic', None, 'mp',

992 ('a', 'b', 'alpha'), [['OBL', 'GYHCH1X', 'GYHCH1XG', None]],

993 ndim=2)

994class OBL(BravaisLattice):

995 def __init__(self, a, b, alpha, **kwargs):

996 check_rect(a, b)

997 if alpha >= 90:

998 raise UnconventionalLattice(

999 f'Expected alpha < 90, got alpha={alpha}')

1000 super().__init__(a=a, b=b, alpha=alpha, **kwargs)

1001

1002 def _cell(self, a, b, alpha):

1003 cosa = np.cos(alpha * _degrees)

1004 sina = np.sin(alpha * _degrees)

1005

1006 return np.array([[a, 0, 0],

1007 [b * cosa, b * sina, 0],

1008 [0., 0., 0.]])

1009

1010 def _special_points(self, a, b, alpha, variant):

1011 cosa = np.cos(alpha * _degrees)

1012 eta = (1 - a * cosa / b) / (2 * np.sin(alpha * _degrees)**2)

1013 nu = .5 - eta * b * cosa / a

1014

1015 points = [[0, 0, 0],

1016 [0, 0.5, 0],

1017 [eta, 1 - nu, 0],

1018 [.5, .5, 0],

1019 [1 - eta, nu, 0],

1020 [.5, 0, 0]]

1021

1022 return points

1023

1024

1025@bravaisclass('primitive hexagonal', 'hexagonal', None, 'hp', 'a',

1026 [['HEX2D', 'GMK', 'GMKG',

1027 get_subset_points('GMK',

1028 sc_special_points['hexagonal'])]],

1029 ndim=2)

1030class HEX2D(BravaisLattice):

1031 def __init__(self, a, **kwargs):

1032 super().__init__(a=a, **kwargs)

1033

1034 def _cell(self, a):

1035 x = 0.5 * np.sqrt(3)

1036 return np.array([[a, 0, 0],

1037 [-0.5 * a, x * a, 0],

1038 [0., 0., 0.]])

1039

1040

1041def check_rect(a, b):

1042 if a >= b:

1043 raise UnconventionalLattice(f'Expected a < b, got a={a}, b={b}')

1044

1045

1046@bravaisclass('primitive rectangular', 'orthorhombic', None, 'op', 'ab',

1047 [['RECT', 'GXSY', 'GXSYGS',

1048 get_subset_points('GXSY',

1049 sc_special_points['orthorhombic'])]],

1050 ndim=2)

1051class RECT(BravaisLattice):

1052 def __init__(self, a, b, **kwargs):

1053 check_rect(a, b)

1054 super().__init__(a=a, b=b, **kwargs)

1055

1056 def _cell(self, a, b):

1057 return np.array([[a, 0, 0],

1058 [0, b, 0],

1059 [0, 0, 0.]])

1060

1061

1062@bravaisclass('centred rectangular', 'orthorhombic', None, 'oc',

1063 ('a', 'alpha'), [['CRECT', 'GXA1Y', 'GXA1YG', None]], ndim=2)

1064class CRECT(BravaisLattice):

1065 def __init__(self, a, alpha, **kwargs):

1066 # It would probably be better to define the CRECT cell

1067 # by (a, b) rather than (a, alpha). Then we can require a < b

1068 # like in ordinary RECT.

1069 #

1070 # In 3D, all lattices in the same family generally take

1071 # identical parameters.

1072 if alpha >= 90:

1073 raise UnconventionalLattice(

1074 f'Expected alpha < 90. Got alpha={alpha}')

1075 super().__init__(a=a, alpha=alpha, **kwargs)

1076

1077 def _cell(self, a, alpha):

1078 x = np.cos(alpha * _degrees)

1079 y = np.sin(alpha * _degrees)

1080 return np.array([[a, 0, 0],

1081 [a * x, a * y, 0],

1082 [0, 0, 0.]])

1083

1084 def _special_points(self, a, alpha, variant):

1085 sina2 = np.sin(alpha / 2 * _degrees)**2

1086 sina = np.sin(alpha * _degrees)**2

1087 eta = sina2 / sina

1088 cosa = np.cos(alpha * _degrees)

1089 xi = eta * cosa

1090

1091 points = [[0, 0, 0],

1092 [eta, - eta, 0],

1093 [0.5 + xi, 0.5 - xi, 0],

1094 [0.5, 0.5, 0]]

1095

1096 return points

1097

1098

1099@bravaisclass('primitive square', 'tetragonal', None, 'tp', ('a',),

1100 [['SQR', 'GMX', 'MGXM',

1101 get_subset_points('GMX', sc_special_points['tetragonal'])]],

1102 ndim=2)

1103class SQR(BravaisLattice):

1104 def __init__(self, a, **kwargs):

1105 super().__init__(a=a, **kwargs)

1106

1107 def _cell(self, a):

1108 return np.array([[a, 0, 0],

1109 [0, a, 0],

1110 [0, 0, 0.]])

1111

1112

1113@bravaisclass('primitive line', 'line', None, '?', ('a',),

1114 [['LINE', 'GX', 'GX', {'G': [0, 0, 0], 'X': [0.5, 0, 0]}]],

1115 ndim=1)

1116class LINE(BravaisLattice):

1117 def __init__(self, a, **kwargs):

1118 super().__init__(a=a, **kwargs)

1119

1120 def _cell(self, a):

1121 return np.array([[a, 0.0, 0.0],

1122 [0.0, 0.0, 0.0],

1123 [0.0, 0.0, 0.0]])

1124

1125

1126def celldiff(cell1, cell2):

1127 """Return a unitless measure of the difference between two cells."""

1128 cell1 = Cell.ascell(cell1).complete()

1129 cell2 = Cell.ascell(cell2).complete()

1130 v1v2 = cell1.volume * cell2.volume

1131 if v1v2 < 1e-10:

1132 # (Proposed cell may be linearly dependent)

1133 return np.inf

1134

1135 scale = v1v2**(-1. / 3.) # --> 1/Ang^2

1136 x1 = cell1 @ cell1.T

1137 x2 = cell2 @ cell2.T

1138 dev = scale * np.abs(x2 - x1).max()

1139 return dev

1140

1141

1142def get_lattice_from_canonical_cell(cell, eps=2e-4):

1143 """Return a Bravais lattice representing the given cell.

1144

1145 This works only for cells that are derived from the standard form

1146 (as generated by lat.tocell()) or rotations thereof.

1147

1148 If the given cell does not resemble the known form of a Bravais

1149 lattice, raise RuntimeError."""

1150 return NormalizedLatticeMatcher(cell, eps).match()

1151

1152

1153class LatticeMatcher:

1154 def __init__(self, cell, pbc=None, *, eps, niggli_eps=None):

1155 self.orig_cell = cell

1156 if pbc is None:

1157 pbc = cell.mask()

1158 self.pbc = cell.any(1) & pbc2pbc(pbc)

1159 self.cell = cell.uncomplete(pbc)

1160 self.eps = eps

1161 if niggli_eps is None:

1162 niggli_eps = eps

1163 self.niggli_cell, self.niggli_op = self.cell.niggli_reduce(

1164 eps=niggli_eps)

1165

1166 # We tabulate the cell's Niggli-mapped versions so we don't need to

1167 # redo any work when the same Niggli-operation appears multiple times

1168 # in the table:

1169 self.memory = {}

1170

1171 def match(self, latname, operations):

1172 matches = []

1173

1174 for op_key in operations:

1175 checker_and_op = self.memory.get(op_key)

1176 if checker_and_op is None:

1177 # op_3x3 is the 3x3 form of the operation that maps

1178 # Niggli cell to AFlow form.

1179 op_3x3 = np.array(op_key).reshape(3, 3)

1180 candidate = Cell(np.linalg.inv(op_3x3.T) @ self.niggli_cell)

1181 checker = NormalizedLatticeMatcher(candidate, eps=self.eps)

1182 self.memory[op_key] = (checker, op_3x3)

1183 else:

1184 checker, op_3x3 = checker_and_op

1185

1186 lat, error = checker.query(latname)

1187 if lat is None or error > self.eps:

1188 continue

1189

1190 # This is the full operation encompassing

1191 # both Niggli reduction of user input and mapping the

1192 # Niggli reduced form to standard (AFlow) form.

1193 op = op_3x3 @ np.linalg.inv(self.niggli_op)

1194 matches.append(Match(lat, op, error))

1195

1196 return matches

1197

1198

1199@dataclass

1200class Match:

1201 lat: BravaisLattice

1202 op: np.ndarray

1203 error: float

1204

1205 @cached_property

1206 def orthogonality_defect(self) -> float:

1207 cell = self.lat.tocell().complete()

1208 return np.prod(cell.lengths()) / cell.volume

1209

1210

1211def identify_lattice(cell, eps=2e-4, *, pbc=True):

1212 """Find Bravais lattice representing this cell.

1213

1214 Returns Bravais lattice object representing the cell along with

1215 an operation that, applied to the cell, yields the same lengths

1216 and angles as the Bravais lattice object."""

1217 from ase.geometry.bravais_type_engine import niggli_op_table

1218

1219 matcher = LatticeMatcher(cell, pbc, eps=eps)

1220

1221 # We loop through the most symmetric kinds (CUB etc.) and return

1222 # the first one we find:

1223 for latname in lattice_check_orders[matcher.cell.rank]:

1224 # There may be multiple Niggli operations that produce valid

1225 # lattices, at least for MCL. In that case we will pick the

1226 # one whose angle is closest to 90, but it means we cannot

1227 # just return the first one we find so we must remember them:

1228 matches = matcher.match(latname, niggli_op_table[latname])

1229

1230 if not matches:

1231 continue # Move to next Bravais lattice

1232

1233 best = min(matches, key=lambda match: match.orthogonality_defect)

1234

1235 if matcher.cell.rank == 2 and best.op[2, 2] < 0:

1236 best.op = flip_2d_handedness(best.op)

1237

1238 return best.lat, best.op

1239

1240 raise RuntimeError('Failed to recognize lattice')

1241

1242

1243def flip_2d_handedness(op):

1244 # The 3x3 operation may flip the z axis, but then the x/y

1245 # components are necessarily also left-handed which

1246 # means a defacto left-handed 2D bandpath.

1247 #

1248 # We repair this by applying an operation that unflips the

1249 # z axis and interchanges x/y:

1250 repair_op = np.array([[0, 1, 0], [1, 0, 0], [0, 0, -1]])

1251 return repair_op @ op

1252

1253

1254# Map of number of dimensions to order in which to check lattices.

1255# We check most symmetric lattices first, so we can terminate early

1256# when a match is found.

1257#

1258# The check order is slightly different than elsewhere listed order,

1259# as we need to check HEX/RHL before the ORCx family.

1260lattice_check_orders = {

1261 1: ['LINE'],

1262 2: ['SQR', 'RECT', 'HEX2D', 'CRECT', 'OBL'],

1263 3: ['CUB', 'FCC', 'BCC', 'TET', 'BCT', 'HEX', 'RHL',

1264 'ORC', 'ORCF', 'ORCI', 'ORCC', 'MCL', 'MCLC', 'TRI']}

1265

1266

1267class NormalizedLatticeMatcher:

1268 # This class checks that a candidate cell matches the normalized

1269 # form of a Bravais lattice. It's used internally by the LatticeMatcher.

1270

1271 def __init__(self, cell, eps=2e-4):

1272 """Generate Bravais lattices that look (or not) like the given cell.

1273

1274 The cell must be reduced to canonical form, i.e., it must

1275 be possible to produce a cell with the same lengths and angles

1276 by directly through one of the Bravais lattice classes.

1277

1278 Generally for internal use (this module).

1279

1280 For each of the 14 Bravais lattices, this object can produce

1281 a lattice object which represents the same cell, or None if

1282 the tolerance eps is not met."""

1283 self.cell = cell

1284 self.eps = eps

1285

1286 self.cellpar = cell.cellpar()

1287 self.lengths = self.A, self.B, self.C = self.cellpar[:3]

1288 self.angles = self.cellpar[3:]

1289

1290 # Use a 'neutral' length for checking cubic lattices

1291 self.A0 = self.lengths.mean()

1292

1293 # Vector of the diagonal and then off-diagonal dot products:

1294 # [a1 · a1, a2 · a2, a3 · a3, a2 · a3, a3 · a1, a1 · a2]

1295 self.prods = (cell @ cell.T).flat[[0, 4, 8, 5, 2, 1]]

1296

1297 def _check(self, latcls, *args):

1298 if any(arg <= 0 for arg in args):

1299 return None

1300

1301 try:

1302 return latcls(*args)

1303 except UnconventionalLattice:

1304 return None

1305

1306 def match(self):

1307 """Match cell against all lattices, returning most symmetric match.

1308

1309 Returns the lattice object. Raises RuntimeError on failure."""

1310 for name in lattice_check_orders[self.cell.rank]:

1311 lat, err = self.query(name)

1312 if lat and err < self.eps:

1313 return lat

1314

1315 raise RuntimeError('Could not find lattice type for cell '

1316 'with lengths and angles {}'

1317 .format(self.cell.cellpar().tolist()))

1318

1319 def query(self, latname):

1320 """Match cell against named Bravais lattice.

1321

1322 Return lattice object on success, None on failure."""

1323 meth = getattr(self, latname)

1324

1325 lat = meth()

1326 if lat is None:

1327 return None, None

1328

1329 newcell = lat.tocell()

1330 err = celldiff(self.cell, newcell)

1331 return lat, err

1332

1333 def LINE(self):

1334 return self._check(LINE, self.lengths[0])

1335

1336 def SQR(self):

1337 return self._check(SQR, self.lengths[0])

1338

1339 def RECT(self):

1340 return self._check(RECT, *self.lengths[:2])

1341

1342 def CRECT(self):

1343 return self._check(CRECT, self.lengths[0], self.angles[2])

1344

1345 def HEX2D(self):

1346 return self._check(HEX2D, self.lengths[0])

1347

1348 def OBL(self):

1349 return self._check(OBL, *self.lengths[:2], self.angles[2])

1350

1351 def CUB(self):

1352 # These methods (CUB, FCC, ...) all return a lattice object if

1353 # it matches, else None.

1354 return self._check(CUB, self.A0)

1355

1356 def FCC(self):

1357 return self._check(FCC, np.sqrt(2) * self.A0)

1358

1359 def BCC(self):

1360 return self._check(BCC, 2.0 * self.A0 / np.sqrt(3))

1361

1362 def TET(self):

1363 return self._check(TET, self.A, self.C)

1364

1365 def _bct_orci_lengths(self):

1366 # Coordinate-system independent relation for BCT and ORCI

1367 # standard cells:

1368 # a1 · a1 + a2 · a3 == a² / 2

1369 # a2 · a2 + a3 · a1 == a² / 2 (BCT)

1370 # == b² / 2 (ORCI)

1371 # a3 · a3 + a1 · a2 == c² / 2

1372 # We use these to get a, b, and c in those cases.

1373 prods = self.prods

1374 lengthsqr = 2.0 * (prods[:3] + prods[3:])

1375 if any(lengthsqr < 0):

1376 return None

1377 return np.sqrt(lengthsqr)

1378

1379 def BCT(self):

1380 lengths = self._bct_orci_lengths()

1381 if lengths is None:

1382 return None

1383 return self._check(BCT, lengths[0], lengths[2])

1384

1385 def HEX(self):

1386 return self._check(HEX, self.A, self.C)

1387

1388 def RHL(self):

1389 return self._check(RHL, self.A, self.angles[0])

1390

1391 def ORC(self):

1392 return self._check(ORC, *self.lengths)

1393

1394 def ORCF(self):

1395 # ORCF standard cell:

1396 # a2 · a3 = a²/4

1397 # a3 · a1 = b²/4

1398 # a1 · a2 = c²/4

1399 prods = self.prods

1400 if all(prods[3:] > 0):

1401 orcf_abc = 2 * np.sqrt(prods[3:])

1402 return self._check(ORCF, *orcf_abc)

1403

1404 def ORCI(self):

1405 lengths = self._bct_orci_lengths()

1406 if lengths is None:

1407 return None

1408 return self._check(ORCI, *lengths)

1409

1410 def _orcc_ab(self):

1411 # ORCC: a1 · a1 + a2 · a3 = a²/2

1412 # a2 · a2 - a2 · a3 = b²/2

1413 prods = self.prods

1414 orcc_sqr_ab = np.empty(2)

1415 orcc_sqr_ab[0] = 2.0 * (prods[0] + prods[5])

1416 orcc_sqr_ab[1] = 2.0 * (prods[1] - prods[5])

1417 if all(orcc_sqr_ab > 0):

1418 return np.sqrt(orcc_sqr_ab)

1419

1420 def ORCC(self):

1421 orcc_lengths_ab = self._orcc_ab()

1422 if orcc_lengths_ab is None:

1423 return None

1424 return self._check(ORCC, *orcc_lengths_ab, self.C)

1425

1426 def MCL(self):

1427 return self._check(MCL, *self.lengths, self.angles[0])

1428

1429 def MCLC(self):

1430 # MCLC is similar to ORCC:

1431 orcc_ab = self._orcc_ab()

1432 if orcc_ab is None:

1433 return None

1434

1435 prods = self.prods

1436 C = self.C

1437 mclc_a, mclc_b = orcc_ab[::-1] # a, b reversed wrt. ORCC

1438 mclc_cosa = 2.0 * prods[3] / (mclc_b * C)

1439 if -1 < mclc_cosa < 1:

1440 mclc_alpha = np.arccos(mclc_cosa) * 180 / np.pi

1441 if mclc_b > C:

1442 # XXX Temporary fix for certain otherwise

1443 # unrecognizable lattices.

1444 #

1445 # This error could happen if the input lattice maps to

1446 # something just outside the domain of conventional

1447 # lattices (less than the tolerance). Our solution is to

1448 # propose a nearby conventional lattice instead, which

1449 # will then be accepted if it's close enough.

1450 mclc_b = 0.5 * (mclc_b + C)

1451 C = mclc_b

1452 return self._check(MCLC, mclc_a, mclc_b, C, mclc_alpha)

1453

1454 def TRI(self):

1455 return self._check(TRI, *self.cellpar)

1456

1457

1458def match_to_lattice(cell, latname: str, pbc=None) -> list[Match]:

1459 """Match cell to a particular Bravais lattice.

1460

1461 For a given input cell, attempt to represent it as a particular

1462 lattice, such as 'FCC', 'ORC', 'ORCI', etc. Return a list of

1463 :class:`ase.lattice.Match` objects which can be compared

1464 so as to get the best match according to metric tensor error

1465 with respect to the input cell or orthogonality defect.

1466 """

1467 from ase.geometry.bravais_type_engine import niggli_op_table

1468

1469 cell = Cell.ascell(cell)

1470 matcher = LatticeMatcher(

1471 cell=cell,

1472 pbc=pbc,

1473 eps=1e6, # We want all candidates and apply eps in our own way

1474 niggli_eps=1e-5,

1475 )

1476

1477 # We should probably sort the matches by "quality".

1478 # Quality must optimize both fitting error and

1479 # orthogonality defect, so some weighting scheme is necessary,

1480 # e.g. minimize one after filtering out the other.

1481 #

1482 # best = min(matches, key=lambda match: match.orthogonality_defect)

1483 #

1484 # Maybe we should take a key callable as argument.

1485 return matcher.match(latname, niggli_op_table[latname])

1486

1487

1488def all_variants():

1489 """For testing and examples; yield all variants of all lattices."""

1490 a, b, c = 3., 4., 5.

1491 alpha = 55.0

1492 yield CUB(a)

1493 yield FCC(a)

1494 yield BCC(a)

1495 yield TET(a, c)

1496 bct1 = BCT(2 * a, c)

1497 bct2 = BCT(a, c)

1498 assert bct1.variant == 'BCT1'

1499 assert bct2.variant == 'BCT2'

1500

1501 yield bct1

1502 yield bct2

1503

1504 yield ORC(a, b, c)

1505

1506 a0 = np.sqrt(1.0 / (1 / b**2 + 1 / c**2))

1507 orcf1 = ORCF(0.5 * a0, b, c)

1508 orcf2 = ORCF(1.2 * a0, b, c)

1509 orcf3 = ORCF(a0, b, c)

1510 assert orcf1.variant == 'ORCF1'

1511 assert orcf2.variant == 'ORCF2'

1512 assert orcf3.variant == 'ORCF3'

1513 yield orcf1

1514 yield orcf2

1515 yield orcf3

1516

1517 yield ORCI(a, b, c)

1518 yield ORCC(a, b, c)

1519

1520 yield HEX(a, c)

1521

1522 rhl1 = RHL(a, alpha=55.0)

1523 assert rhl1.variant == 'RHL1'

1524 yield rhl1

1525

1526 rhl2 = RHL(a, alpha=105.0)

1527 assert rhl2.variant == 'RHL2'

1528 yield rhl2

1529

1530 # With these lengths, alpha < 65 (or so) would result in a lattice that

1531 # could also be represented with alpha > 65, which is more conventional.

1532 yield MCL(a, b, c, alpha=70.0)

1533

1534 mclc1 = MCLC(a, b, c, 80)

1535 assert mclc1.variant == 'MCLC1'

1536 yield mclc1

1537 # mclc2 has same special points as mclc1

1538

1539 mclc3 = MCLC(1.8 * a, b, c * 2, 80)

1540 assert mclc3.variant == 'MCLC3'

1541 yield mclc3

1542 # mclc4 has same special points as mclc3

1543

1544 # XXX We should add MCLC2 and MCLC4 as well.

1545

1546 mclc5 = MCLC(b, b, 1.1 * b, 70)

1547 assert mclc5.variant == 'MCLC5'

1548 yield mclc5

1549

1550 def get_tri(kcellpar):

1551 # We build the TRI lattices from cellpars of reciprocal cell

1552 icell = Cell.fromcellpar(kcellpar)

1553 cellpar = Cell(4 * icell.reciprocal()).cellpar()

1554 return TRI(*cellpar)

1555

1556 tri1a = get_tri([1., 1.2, 1.4, 120., 110., 100.])

1557 assert tri1a.variant == 'TRI1a'

1558 yield tri1a

1559

1560 tri1b = get_tri([1., 1.2, 1.4, 50., 60., 70.])

1561 assert tri1b.variant == 'TRI1b'

1562 yield tri1b

1563

1564 tri2a = get_tri([1., 1.2, 1.4, 120., 110., 90.])

1565 assert tri2a.variant == 'TRI2a'

1566 yield tri2a

1567 tri2b = get_tri([1., 1.2, 1.4, 50., 60., 90.])

1568 assert tri2b.variant == 'TRI2b'

1569 yield tri2b

1570

1571 # Choose an OBL lattice that round-trip-converts to itself.

1572 # The default a/b/alpha parameters result in another representation

1573 # of the same lattice.

1574 yield OBL(a=3.0, b=3.35, alpha=77.85)

1575 yield RECT(a, b)

1576 yield CRECT(a, alpha=alpha)

1577 yield HEX2D(a)

1578 yield SQR(a)

1579 yield LINE(a)

Coverage for ase / lattice / init.py: 97.13%

767 statements