Coverage for ase/lattice/__init_

1# fmt: off

3from abc import ABC, abstractmethod

4from dataclasses import dataclass

5from functools import cached_property

7import numpy as np

9from ase.cell import Cell

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

11from ase.utils import pbc2pbc

13_degrees = np.pi / 180

16def lattice_attr(name):

17 def getter(self):

18 try:

19 return self._parameters[name]

20 except KeyError:

21 raise AttributeError(name) from None

23 getter.__name__ = name

24 return property(getter)

27class BravaisLattice(ABC):

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

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

31 five 2D classes.

33 Each class stores basic static information:

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

36 >>> FCC.name

37 'FCC'

38 >>> FCC.longname

39 'face-centred cubic'

40 >>> FCC.pearson_symbol

41 'cF'

42 >>> MCL.parameters

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

45 Each class can be instantiated with the specific lattice parameters

46 that apply to that lattice:

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

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

51 """

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

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

54 #

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

56 # not None, unfortunately.

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

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

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

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

61 # 'BCT2': <variant object>}

62 ndim: int = 0

63 conventional_cls: str | None = None

65 def __init__(self, **kwargs):

66 p = {}

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

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

69 p[k] = float(v)

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

71 self._parameters = p

72 self._eps = eps

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

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

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

77 else:

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

79 self._variant = self.variants[name]

81 @property

82 def variant(self) -> str:

83 """Return name of lattice variant.

85 >>> from ase.lattice import BCT

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

88 'BCT2'

89 """

90 return self._variant.name

92 def vars(self) -> dict[str, float]:

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

94 return dict(self._parameters)

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

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

98 cls = bravais_lattices[self.conventional_cls]

99 return cls(**self._parameters)

100

101 def tocell(self) -> Cell:

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

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

104 return Cell(cell)

105

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

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

108

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

110 # (Just a brute-force implementation)

111 cell = self.tocell()

112 return cell.cellpar()

113

114 @property

115 def special_path(self) -> str:

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

117

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

119 'GXYSGZS1NPY1Z,XP'

120 """

121 return self._variant.special_path

122

123 @property

124 def special_point_names(self) -> list[str]:

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

126

127 >>> from ase.lattice import BCT

128

129 >>> BCT(3, 5).special_point_names

130 ['G', 'N', 'P', 'S', 'S1', 'X', 'Y', 'Y1', 'Z']

131 """

132 labels = parse_path_string(self._variant.special_point_names)

133 assert len(labels) == 1 # list of lists

134 return labels[0]

135

136 def get_special_points_array(self) -> np.ndarray:

137 """Return all special points for this lattice as an array.

138

139 Ordering is consistent with special_point_names."""

140 if self._variant.special_points is not None:

141 # Fixed dictionary of special points

142 d = self.get_special_points()

143 labels = self.special_point_names

144 assert len(d) == len(labels)

145 points = np.empty((len(d), 3))

146 for i, label in enumerate(labels):

147 points[i] = d[label]

148 return points

149

150 # Special points depend on lattice parameters:

151 points = self._special_points(variant=self._variant,

152 **self._parameters)

153 assert len(points) == len(self.special_point_names)

154 return np.array(points)

155

156 def get_special_points(self) -> dict[str, np.ndarray]:

157 """Return a dictionary of named special k-points for this lattice."""

158 if self._variant.special_points is not None:

159 return self._variant.special_points

160

161 labels = self.special_point_names

162 points = self.get_special_points_array()

163

164 return dict(zip(labels, points))

165

166 def plot_bz(self, path=None, special_points=None, **plotkwargs):

167 """Plot the reciprocal cell and default bandpath."""

168 # Create a generic bandpath (no interpolated kpoints):

169 bandpath = self.bandpath(path=path, special_points=special_points,

170 npoints=0)

171 return bandpath.plot(dimension=self.ndim, **plotkwargs)

172

173 def bandpath(self, path=None, npoints=None, special_points=None,

174 density=None) -> BandPath:

175 """Return a :class:`~ase.dft.kpoints.BandPath` for this lattice.

176

177 See :meth:`ase.cell.Cell.bandpath` for description of parameters.

178

179 >>> from ase.lattice import BCT

180

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

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

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

184

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

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

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

188 the kpoints may not correspond to your particular cell.

189

190 """

191 if special_points is None:

192 special_points = self.get_special_points()

193

194 if path is None:

195 path = self._variant.special_path

196 elif not isinstance(path, str):

197 from ase.dft.kpoints import resolve_custom_points

198 path, special_points = resolve_custom_points(path,

199 special_points,

200 self._eps)

201

202 cell = self.tocell()

203 bandpath = BandPath(cell=cell, path=path,

204 special_points=special_points)

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

206

207 @abstractmethod

208 def _cell(self, **kwargs):

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

210

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

212

213 def _special_points(self, **kwargs):

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

215

216 Subclasses typically return a nested list.

217

218 Ordering must be same as kpoint labels.

219

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

221 raise NotImplementedError

222

223 def _variant_name(self, **kwargs):

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

225

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

227 raise NotImplementedError

228

229 def __format__(self, spec):

230 tokens = []

231 if not spec:

232 spec = '.6g'

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

234

235 for name in self.parameters:

236 value = self._parameters[name]

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

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

239

240 def __str__(self) -> str:

241 return self.__format__('')

242

243 def __repr__(self) -> str:

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

245

246 def description(self) -> str:

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

248 points = self.get_special_points()

249 labels = self.special_point_names

250

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

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

253 for label in labels])

254

255 string = """\

256{repr}

257 {variant}

258 Special point coordinates:

259{special_points}

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

261 variant=self._variant,

262 special_points=coordstring)

263 return string

264

265 @classmethod

266 def type_description(cls):

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

268 desc = """\

269Lattice name: {name}

270 Long name: {longname}

271 Parameters: {parameters}

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

273

274 chunks = [desc]

275 for name in cls.variant_names:

276 var = cls.variants[name]

277 txt = str(var)

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

279 lines.append('')

280 chunks.extend(lines)

281

282 return '\n'.join(chunks)

283

284

285class Variant:

286 variant_desc = """\

287Variant name: {name}

288 Special point names: {special_point_names}

289 Default path: {special_path}

290"""

291

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

293 special_points=None):

294 self.name = name

295 self.special_point_names = special_point_names

296 self.special_path = special_path

297 if special_points is not None:

298 special_points = dict(special_points)

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

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

301 self.special_points = special_points

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

303 # (easier to transform that way)

304

305 def __str__(self) -> str:

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

307

308

309bravais_names = []

310bravais_lattices = {}

311bravais_classes = {}

312

313

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

315 parameters, variants, ndim=3):

316 """Decorator for Bravais lattice classes.

317

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

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

320

321 def decorate(cls):

322 btype = cls.__name__

323 cls.name = btype

324 cls.longname = longname

325 cls.crystal_family = crystal_family

326 cls.lattice_system = lattice_system

327 cls.pearson_symbol = pearson_symbol

328 cls.parameters = tuple(parameters)

329 cls.variant_names = []

330 cls.variants = {}

331 cls.ndim = ndim

332

333 for parameter in parameters:

334 setattr(cls, parameter, lattice_attr(parameter))

335

336 for [name, special_point_names, special_path,

337 special_points] in variants:

338 cls.variant_names.append(name)

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

340 special_path, special_points)

341

342 # Register in global list and dictionary

343 bravais_names.append(btype)

344 bravais_lattices[btype] = cls

345 bravais_classes[pearson_symbol] = cls

346 return cls

347

348 return decorate

349

350

351# Common mappings of primitive to conventional cells:

352_identity = np.identity(3, int)

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

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

355

356

357class UnconventionalLattice(ValueError):

358 pass

359

360

361class Cubic(BravaisLattice):

362 """Abstract class for cubic lattices."""

363 conventional_cls = 'CUB'

364

365 def __init__(self, a):

366 super().__init__(a=a)

367

368

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

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

371class CUB(Cubic):

372 conventional_cellmap = _identity

373

374 def _cell(self, a):

375 return a * np.eye(3)

376

377

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

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

380class FCC(Cubic):

381 conventional_cellmap = _bcc_map

382

383 def _cell(self, a):

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

385

386

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

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

389class BCC(Cubic):

390 conventional_cellmap = _fcc_map

391

392 def _cell(self, a):

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

394

395

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

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

398 sc_special_points['tetragonal']]])

399class TET(BravaisLattice):

400 conventional_cls = 'TET'

401 conventional_cellmap = _identity

402

403 def __init__(self, a, c):

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

405

406 def _cell(self, a, c):

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

408

409

410# XXX in BCT2 we use S for Sigma.

411# Also in other places I think

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

413 'ac',

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

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

416class BCT(BravaisLattice):

417 conventional_cls = 'TET'

418 conventional_cellmap = _fcc_map

419

420 def __init__(self, a, c):

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

422

423 def _cell(self, a, c):

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

425

426 def _variant_name(self, a, c):

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

428

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

430 a2 = a * a

431 c2 = c * c

432

433 assert variant.name in self.variants

434

435 if variant.name == 'BCT1':

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

437 points = [[0, 0, 0],

438 [-.5, .5, .5],

439 [0., .5, 0.],

440 [.25, .25, .25],

441 [0., 0., .5],

442 [eta, eta, -eta],

443 [-eta, 1 - eta, eta]]

444 else:

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

446 zeta = 0.5 * a2 / c2

447 points = [[0., .0, 0.],

448 [0., .5, 0.],

449 [.25, .25, .25],

450 [-eta, eta, eta],

451 [eta, 1 - eta, -eta],

452 [0., 0., .5],

453 [-zeta, zeta, .5],

454 [.5, .5, -zeta],

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

456 return points

457

458

459def check_orc(a, b, c):

460 if not a < b < c:

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

462 .format(a, b, c))

463

464

465class Orthorhombic(BravaisLattice):

466 """Abstract class for orthorhombic types."""

467

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

469 check_orc(a, b, c)

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

471

472

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

474 'abc',

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

476 sc_special_points['orthorhombic']]])

477class ORC(Orthorhombic):

478 conventional_cls = 'ORC'

479 conventional_cellmap = _identity

480

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

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

483

484

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

486 'oF', 'abc',

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

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

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

490class ORCF(Orthorhombic):

491 conventional_cls = 'ORC'

492 conventional_cellmap = _bcc_map

493

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

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

496

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

498 a2 = a * a

499 b2 = b * b

500 c2 = c * c

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

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

503

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

505 zeta = xminus

506 eta = xplus

507

508 points = [[0, 0, 0],

509 [.5, .5 + zeta, zeta],

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

511 [.5, .5, .5],

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

513 [0., eta, eta],

514 [1., 1 - eta, 1 - eta],

515 [.5, 0, .5],

516 [.5, .5, 0]]

517 else:

518 assert variant.name == 'ORCF2'

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

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

521 eta = xminus

522

523 points = [[0, 0, 0],

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

525 [.5, .5 + eta, eta],

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

527 [.5 + delta, .5, delta],

528 [.5, .5, .5],

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

530 [phi, .5 + phi, .5],

531 [0., .5, .5],

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

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

534

535 return points

536

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

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

539 if abs(diff) < self._eps:

540 return 'ORCF3'

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

542

543

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

545 'oI', 'abc',

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

547class ORCI(Orthorhombic):

548 conventional_cls = 'ORC'

549 conventional_cellmap = _fcc_map

550

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

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

553

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

555 a2 = a**2

556 b2 = b**2

557 c2 = c**2

558

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

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

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

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

563

564 points = [[0., 0., 0.],

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

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

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

568 [0, .5, 0],

569 [.5, 0, 0],

570 [0., 0., .5],

571 [.25, .25, .25],

572 [-zeta, zeta, zeta],

573 [zeta, 1 - zeta, -zeta],

574 [eta, -eta, eta],

575 [1 - eta, eta, -eta],

576 [.5, .5, -.5]]

577 return points

578

579

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

581 'oC', 'abc',

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

583class ORCC(BravaisLattice):

584 conventional_cls = 'ORC'

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

586

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

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

589 if a >= b:

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

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

592

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

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

595 [0, 0, c]])

596

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

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

599 points = [[0, 0, 0],

600 [zeta, zeta, .5],

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

602 [0, .5, .5],

603 [0, .5, 0],

604 [-.5, .5, .5],

605 [zeta, zeta, 0],

606 [-zeta, 1 - zeta, 0],

607 [-.5, .5, 0],

608 [0, 0, .5]]

609 return points

610

611

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

613 'ac',

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

615 sc_special_points['hexagonal']]])

616class HEX(BravaisLattice):

617 conventional_cls = 'HEX'

618 conventional_cellmap = _identity

619

620 def __init__(self, a, c):

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

622

623 def _cell(self, a, c):

624 x = 0.5 * np.sqrt(3)

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

626 [0., 0., c]])

627

628

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

630 ('a', 'alpha'),

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

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

633class RHL(BravaisLattice):

634 conventional_cls = 'RHL'

635 conventional_cellmap = _identity

636

637 def __init__(self, a, alpha):

638 if alpha >= 120:

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

640 .format(alpha))

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

642

643 def _cell(self, a, alpha):

644 alpha *= np.pi / 180

645 acosa = a * np.cos(alpha)

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

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

648 acosfrac = acosa / acosa2

649 xx = (1 - acosfrac**2)

650 assert xx > 0.0

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

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

653

654 def _variant_name(self, a, alpha):

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

656

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

658 if variant.name == 'RHL1':

659 cosa = np.cos(alpha * _degrees)

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

661 nu = .75 - 0.5 * eta

662 points = [[0, 0, 0],

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

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

665 [.5, .5, 0],

666 [.5, 0, 0],

667 [0, 0, -.5],

668 [eta, nu, nu],

669 [1 - nu, 1 - nu, 1 - eta],

670 [nu, nu, eta - 1],

671 [1 - nu, nu, 0],

672 [nu, 0, -nu],

673 [.5, .5, .5]]

674 else:

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

676 nu = .75 - 0.5 * eta

677 points = [[0, 0, 0],

678 [.5, -.5, 0],

679 [.5, 0, 0],

680 [1 - nu, -nu, 1 - nu],

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

682 [eta, eta, eta],

683 [1 - eta, -eta, -eta],

684 [.5, -.5, .5]]

685 return points

686

687

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

689 if b > c or alpha >= 90:

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

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

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

693

694

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

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

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

698class MCL(BravaisLattice):

699 conventional_cls = 'MCL'

700 conventional_cellmap = _identity

701

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

703 check_mcl(a, b, c, alpha)

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

705

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

707 alpha *= _degrees

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

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

710

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

712 cosa = np.cos(alpha * _degrees)

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

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

715

716 points = [[0, 0, 0],

717 [.5, .5, 0],

718 [0, .5, .5],

719 [.5, 0, .5],

720 [.5, 0, -.5],

721 [.5, .5, .5],

722 [0, eta, 1 - nu],

723 [0, 1 - eta, nu],

724 [0, eta, -nu],

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

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

727 [.5, eta, -nu],

728 [0, .5, 0],

729 [0, 0, .5],

730 [0, 0, -.5],

731 [.5, 0, 0]]

732 return points

733

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

735 check_mcl(a, b, c, alpha)

736 return 'MCL'

737

738

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

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

741 [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

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

743 ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

744 'GYFLI,I1ZF1,NGM', None],

745 ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

746 'GYFHZIF1,H1Y1XGN,MG', None],

747 ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

748 'GYFHZI,H1Y1XGN,MG', None],

749 ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z',

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

751class MCLC(BravaisLattice):

752 conventional_cls = 'MCL'

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

754

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

756 check_mcl(a, b, c, alpha)

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

758

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

760 alpha *= np.pi / 180

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

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

763

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

765 # from ase.geometry.cell import mclc

766 # okay, this is a bit hacky

767

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

769 # Right now we just repeat the code:

770 check_mcl(a, b, c, alpha)

771

772 a2 = a * a

773 b2 = b * b

774 cosa = np.cos(alpha * _degrees)

775 sina = np.sin(alpha * _degrees)

776 sina2 = sina**2

777

778 cell = self.tocell()

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

780

781 kgamma = lengths_angles[-1]

782

783 eps = self._eps

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

785 # the same precision.

786 if abs(kgamma - 90) < eps:

787 variant = 2

788 elif kgamma > 90:

789 variant = 1

790 elif kgamma < 90:

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

792 if abs(num - 1) < eps:

793 variant = 4

794 elif num < 1:

795 variant = 3

796 else:

797 variant = 5

798 variant = 'MCLC' + str(variant)

799 return variant

800

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

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

803

804 a2 = a * a

805 b2 = b * b

806 # c2 = c * c

807 cosa = np.cos(alpha * _degrees)

808 sina = np.sin(alpha * _degrees)

809 sina2 = sina**2

810

811 if variant == 1 or variant == 2:

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

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

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

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

816

817 points = [[0, 0, 0],

818 [.5, 0, 0],

819 [0, -.5, 0],

820 [1 - zeta, 1 - zeta, 1 - eta],

821 [zeta, zeta, eta],

822 [-zeta, -zeta, 1 - eta],

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

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

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

826 [.5, .5, .5],

827 [.5, 0, .5],

828 [1 - psi, psi - 1, 0],

829 [psi, 1 - psi, 0],

830 [psi - 1, -psi, 0],

831 [.5, .5, 0],

832 [-.5, -.5, 0],

833 [0, 0, .5]]

834 elif variant == 3 or variant == 4:

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

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

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

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

839 phi = 1 + zeta - 2 * mu

840 psi = eta - 2 * delta

841

842 points = [[0, 0, 0],

843 [1 - phi, 1 - phi, 1 - psi],

844 [phi, phi - 1, psi],

845 [1 - phi, -phi, 1 - psi],

846 [zeta, zeta, eta],

847 [1 - zeta, -zeta, 1 - eta],

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

849 [.5, -.5, .5],

850 [.5, 0, .5],

851 [.5, 0, 0],

852 [0, -.5, 0],

853 [.5, -.5, 0],

854 [mu, mu, delta],

855 [1 - mu, -mu, -delta],

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

857 [mu, mu - 1, delta],

858 [0, 0, .5]]

859 elif variant == 5:

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

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

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

863 nu = 2 * mu - zeta

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

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

866 rho = 1 - zeta * a2 / b2

867

868 points = [[0, 0, 0],

869 [nu, nu, omega],

870 [1 - nu, 1 - nu, 1 - omega],

871 [nu, nu - 1, omega],

872 [zeta, zeta, eta],

873 [1 - zeta, -zeta, 1 - eta],

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

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

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

877 [.5, .5, .5],

878 [.5, 0, .5],

879 [.5, 0, 0],

880 [0, -.5, 0],

881 [.5, -.5, 0],

882 [mu, mu, delta],

883 [1 - mu, -mu, -delta],

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

885 [mu, mu - 1, delta],

886 [0, 0, .5]]

887

888 return points

889

890

891tri_angles_explanation = """\

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

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

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

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

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

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

898

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

900

901

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

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

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

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

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

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

908class TRI(BravaisLattice):

909 conventional_cls = 'TRI'

910 conventional_cellmap = _identity

911

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

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

914 gamma=gamma)

915

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

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

918 singamma = np.sin(gamma * _degrees)

919 cosgamma = np.cos(gamma * _degrees)

920 cosbeta = np.cos(beta * _degrees)

921 cosalpha = np.cos(alpha * _degrees)

922 a3x = c * cosbeta

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

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

925 + 2 * cosalpha * cosbeta * cosgamma)

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

927 [a3x, a3y, a3z]])

928

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

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

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

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

933

934 def raise_unconventional():

935 raise UnconventionalLattice(tri_angles_explanation

936 .format(*kangles))

937

938 eps = self._eps

939 if abs(kgamma - 90) < eps:

940 if kalpha > 90 and kbeta > 90:

941 var = '2a'

942 elif kalpha < 90 and kbeta < 90:

943 var = '2b'

944 else:

945 # Is this possible? Maybe due to epsilon

946 raise_unconventional()

947 elif all(kangles > 90):

948 if kgamma > min(kangles):

949 raise_unconventional()

950 var = '1a'

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

952 if kgamma < max(kangles):

953 raise_unconventional()

954 var = '1b'

955 else:

956 raise_unconventional()

957

958 return 'TRI' + var

959

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

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

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

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

964 points = [[0., 0., 0.],

965 [.5, .5, 0],

966 [0, .5, .5],

967 [.5, 0, .5],

968 [.5, .5, .5],

969 [.5, 0, 0],

970 [0, .5, 0],

971 [0, 0, .5]]

972 else:

973 points = [[0, 0, 0],

974 [.5, -.5, 0],

975 [0, 0, .5],

976 [-.5, -.5, .5],

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

978 [0, -0.5, 0],

979 [.5, 0, 0],

980 [-.5, 0, .5]]

981

982 return points

983

984

985def get_subset_points(names, points):

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

987 return newpoints

988

989

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

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

992 ndim=2)

993class OBL(BravaisLattice):

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

995 check_rect(a, b)

996 if alpha >= 90:

997 raise UnconventionalLattice(

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

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

1000

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

1002 cosa = np.cos(alpha * _degrees)

1003 sina = np.sin(alpha * _degrees)

1004

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

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

1007 [0., 0., 0.]])

1008

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

1010 cosa = np.cos(alpha * _degrees)

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

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

1013

1014 points = [[0, 0, 0],

1015 [0, 0.5, 0],

1016 [eta, 1 - nu, 0],

1017 [.5, .5, 0],

1018 [1 - eta, nu, 0],

1019 [.5, 0, 0]]

1020

1021 return points

1022

1023

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

1025 [['HEX2D', 'GMK', 'GMKG',

1026 get_subset_points('GMK',

1027 sc_special_points['hexagonal'])]],

1028 ndim=2)

1029class HEX2D(BravaisLattice):

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

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

1032

1033 def _cell(self, a):

1034 x = 0.5 * np.sqrt(3)

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

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

1037 [0., 0., 0.]])

1038

1039

1040def check_rect(a, b):

1041 if a >= b:

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

1043

1044

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

1046 [['RECT', 'GXSY', 'GXSYGS',

1047 get_subset_points('GXSY',

1048 sc_special_points['orthorhombic'])]],

1049 ndim=2)

1050class RECT(BravaisLattice):

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

1052 check_rect(a, b)

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

1054

1055 def _cell(self, a, b):

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

1057 [0, b, 0],

1058 [0, 0, 0.]])

1059

1060

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

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

1063class CRECT(BravaisLattice):

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

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

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

1067 # like in ordinary RECT.

1068 #

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

1070 # identical parameters.

1071 if alpha >= 90:

1072 raise UnconventionalLattice(

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

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

1075

1076 def _cell(self, a, alpha):

1077 x = np.cos(alpha * _degrees)

1078 y = np.sin(alpha * _degrees)

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

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

1081 [0, 0, 0.]])

1082

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

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

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

1086 eta = sina2 / sina

1087 cosa = np.cos(alpha * _degrees)

1088 xi = eta * cosa

1089

1090 points = [[0, 0, 0],

1091 [eta, - eta, 0],

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

1093 [0.5, 0.5, 0]]

1094

1095 return points

1096

1097

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

1099 [['SQR', 'GMX', 'MGXM',

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

1101 ndim=2)

1102class SQR(BravaisLattice):

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

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

1105

1106 def _cell(self, a):

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

1108 [0, a, 0],

1109 [0, 0, 0.]])

1110

1111

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

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

1114 ndim=1)

1115class LINE(BravaisLattice):

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

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

1118

1119 def _cell(self, a):

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

1121 [0.0, 0.0, 0.0],

1122 [0.0, 0.0, 0.0]])

1123

1124

1125def celldiff(cell1, cell2):

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

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

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

1129 v1v2 = cell1.volume * cell2.volume

1130 if v1v2 < 1e-10:

1131 # (Proposed cell may be linearly dependent)

1132 return np.inf

1133

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

1135 x1 = cell1 @ cell1.T

1136 x2 = cell2 @ cell2.T

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

1138 return dev

1139

1140

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

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

1143

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

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

1146

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

1148 lattice, raise RuntimeError."""

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

1150

1151

1152class LatticeMatcher:

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

1154 self.orig_cell = cell

1155 if pbc is None:

1156 pbc = cell.mask()

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

1158 self.cell = cell.uncomplete(pbc)

1159 self.eps = eps

1160 if niggli_eps is None:

1161 niggli_eps = eps

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

1163 eps=niggli_eps)

1164

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

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

1167 # in the table:

1168 self.memory = {}

1169

1170 def match(self, latname, operations):

1171 matches = []

1172

1173 for op_key in operations:

1174 checker_and_op = self.memory.get(op_key)

1175 if checker_and_op is None:

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

1177 # Niggli cell to AFlow form.

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

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

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

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

1182 else:

1183 checker, op_3x3 = checker_and_op

1184

1185 lat, error = checker.query(latname)

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

1187 continue

1188

1189 # This is the full operation encompassing

1190 # both Niggli reduction of user input and mapping the

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

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

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

1194

1195 return matches

1196

1197

1198@dataclass

1199class Match:

1200 lat: BravaisLattice

1201 op: np.ndarray

1202 error: float

1203

1204 @cached_property

1205 def orthogonality_defect(self) -> float:

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

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

1208

1209

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

1211 """Find Bravais lattice representing this cell.

1212

1213 Returns Bravais lattice object representing the cell along with

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

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

1216 from ase.geometry.bravais_type_engine import niggli_op_table

1217

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

1219

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

1221 # the first one we find:

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

1223 # There may be multiple Niggli operations that produce valid

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

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

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

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

1228

1229 if not matches:

1230 continue # Move to next Bravais lattice

1231

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

1233

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

1235 best.op = flip_2d_handedness(best.op)

1236

1237 return best.lat, best.op

1238

1239 raise RuntimeError('Failed to recognize lattice')

1240

1241

1242def flip_2d_handedness(op):

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

1244 # components are necessarily also left-handed which

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

1246 #

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

1248 # z axis and interchanges x/y:

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

1250 return repair_op @ op

1251

1252

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

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

1255# when a match is found.

1256#

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

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

1259lattice_check_orders = {

1260 1: ['LINE'],

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

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

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

1264

1265

1266class NormalizedLatticeMatcher:

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

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

1269

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

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

1272

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

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

1275 by directly through one of the Bravais lattice classes.

1276

1277 Generally for internal use (this module).

1278

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

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

1281 the tolerance eps is not met."""

1282 self.cell = cell

1283 self.eps = eps

1284

1285 self.cellpar = cell.cellpar()

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

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

1288

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

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

1291

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

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

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

1295

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

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

1298 return None

1299

1300 try:

1301 return latcls(*args)

1302 except UnconventionalLattice:

1303 return None

1304

1305 def match(self):

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

1307

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

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

1310 lat, err = self.query(name)

1311 if lat and err < self.eps:

1312 return lat

1313

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

1315 'with lengths and angles {}'

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

1317

1318 def query(self, latname):

1319 """Match cell against named Bravais lattice.

1320

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

1322 meth = getattr(self, latname)

1323

1324 lat = meth()

1325 if lat is None:

1326 return None, None

1327

1328 newcell = lat.tocell()

1329 err = celldiff(self.cell, newcell)

1330 return lat, err

1331

1332 def LINE(self):

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

1334

1335 def SQR(self):

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

1337

1338 def RECT(self):

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

1340

1341 def CRECT(self):

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

1343

1344 def HEX2D(self):

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

1346

1347 def OBL(self):

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

1349

1350 def CUB(self):

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

1352 # it matches, else None.

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

1354

1355 def FCC(self):

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

1357

1358 def BCC(self):

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

1360

1361 def TET(self):

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

1363

1364 def _bct_orci_lengths(self):

1365 # Coordinate-system independent relation for BCT and ORCI

1366 # standard cells:

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

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

1369 # == b² / 2 (ORCI)

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

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

1372 prods = self.prods

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

1374 if any(lengthsqr < 0):

1375 return None

1376 return np.sqrt(lengthsqr)

1377

1378 def BCT(self):

1379 lengths = self._bct_orci_lengths()

1380 if lengths is None:

1381 return None

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

1383

1384 def HEX(self):

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

1386

1387 def RHL(self):

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

1389

1390 def ORC(self):

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

1392

1393 def ORCF(self):

1394 # ORCF standard cell:

1395 # a2 · a3 = a²/4

1396 # a3 · a1 = b²/4

1397 # a1 · a2 = c²/4

1398 prods = self.prods

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

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

1401 return self._check(ORCF, *orcf_abc)

1402

1403 def ORCI(self):

1404 lengths = self._bct_orci_lengths()

1405 if lengths is None:

1406 return None

1407 return self._check(ORCI, *lengths)

1408

1409 def _orcc_ab(self):

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

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

1412 prods = self.prods

1413 orcc_sqr_ab = np.empty(2)

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

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

1416 if all(orcc_sqr_ab > 0):

1417 return np.sqrt(orcc_sqr_ab)

1418

1419 def ORCC(self):

1420 orcc_lengths_ab = self._orcc_ab()

1421 if orcc_lengths_ab is None:

1422 return None

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

1424

1425 def MCL(self):

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

1427

1428 def MCLC(self):

1429 # MCLC is similar to ORCC:

1430 orcc_ab = self._orcc_ab()

1431 if orcc_ab is None:

1432 return None

1433

1434 prods = self.prods

1435 C = self.C

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

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

1438 if -1 < mclc_cosa < 1:

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

1440 if mclc_b > C:

1441 # XXX Temporary fix for certain otherwise

1442 # unrecognizable lattices.

1443 #

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

1445 # something just outside the domain of conventional

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

1447 # propose a nearby conventional lattice instead, which

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

1449 mclc_b = 0.5 * (mclc_b + C)

1450 C = mclc_b

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

1452

1453 def TRI(self):

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

1455

1456

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

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

1459

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

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

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

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

1464 with respect to the input cell or orthogonality defect.

1465 """

1466 from ase.geometry.bravais_type_engine import niggli_op_table

1467

1468 cell = Cell.ascell(cell)

1469 matcher = LatticeMatcher(

1470 cell=cell,

1471 pbc=pbc,

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

1473 niggli_eps=1e-5,

1474 )

1475

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

1477 # Quality must optimize both fitting error and

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

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

1480 #

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

1482 #

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

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

1485

1486

1487def all_variants():

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

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

1490 alpha = 55.0

1491 yield CUB(a)

1492 yield FCC(a)

1493 yield BCC(a)

1494 yield TET(a, c)

1495 bct1 = BCT(2 * a, c)

1496 bct2 = BCT(a, c)

1497 assert bct1.variant == 'BCT1'

1498 assert bct2.variant == 'BCT2'

1499

1500 yield bct1

1501 yield bct2

1502

1503 yield ORC(a, b, c)

1504

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

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

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

1508 orcf3 = ORCF(a0, b, c)

1509 assert orcf1.variant == 'ORCF1'

1510 assert orcf2.variant == 'ORCF2'

1511 assert orcf3.variant == 'ORCF3'

1512 yield orcf1

1513 yield orcf2

1514 yield orcf3

1515

1516 yield ORCI(a, b, c)

1517 yield ORCC(a, b, c)

1518

1519 yield HEX(a, c)

1520

1521 rhl1 = RHL(a, alpha=55.0)

1522 assert rhl1.variant == 'RHL1'

1523 yield rhl1

1524

1525 rhl2 = RHL(a, alpha=105.0)

1526 assert rhl2.variant == 'RHL2'

1527 yield rhl2

1528

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

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

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

1532

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

1534 assert mclc1.variant == 'MCLC1'

1535 yield mclc1

1536 # mclc2 has same special points as mclc1

1537

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

1539 assert mclc3.variant == 'MCLC3'

1540 yield mclc3

1541 # mclc4 has same special points as mclc3

1542

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

1544

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

1546 assert mclc5.variant == 'MCLC5'

1547 yield mclc5

1548

1549 def get_tri(kcellpar):

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

1551 icell = Cell.fromcellpar(kcellpar)

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

1553 return TRI(*cellpar)

1554

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

1556 assert tri1a.variant == 'TRI1a'

1557 yield tri1a

1558

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

1560 assert tri1b.variant == 'TRI1b'

1561 yield tri1b

1562

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

1564 assert tri2a.variant == 'TRI2a'

1565 yield tri2a

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

1567 assert tri2b.variant == 'TRI2b'

1568 yield tri2b

1569

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

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

1572 # of the same lattice.

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

1574 yield RECT(a, b)

1575 yield CRECT(a, alpha=alpha)

1576 yield HEX2D(a)

1577 yield SQR(a)

1578 yield LINE(a)

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

766 statements