Coverage for /builds/ase/ase/ase/lattice/__init_

1# fmt: off

3from abc import ABC, abstractmethod

4from functools import cached_property

5from typing import Dict, List

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

16class BravaisLattice(ABC):

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

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

20 five 2D classes.

22 Each class stores basic static information:

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

25 >>> FCC.name

26 'FCC'

27 >>> FCC.longname

28 'face-centred cubic'

29 >>> FCC.pearson_symbol

30 'cF'

31 >>> MCL.parameters

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

34 Each class can be instantiated with the specific lattice parameters

35 that apply to that lattice:

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

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

40 """

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

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

43 name = None # e.g. 'CUB', 'BCT', 'ORCF', ...

44 longname = None # e.g. 'cubic', 'body-centred tetragonal', ...

45 parameters = None # e.g. ('a', 'c')

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

47 # 'BCT2': <variant object>}

48 ndim = None

50 def __init__(self, **kwargs):

51 p = {}

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

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

54 p[k] = float(v)

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

56 self._parameters = p

57 self._eps = eps

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

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

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

62 else:

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

64 self._variant = self.variants[name]

66 @property

67 def variant(self) -> str:

68 """Return name of lattice variant.

70 >>> from ase.lattice import BCT

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

73 'BCT2'

74 """

75 return self._variant.name

77 def __getattr__(self, name: str):

78 if name in self._parameters:

79 return self._parameters[name]

80 return self.__getattribute__(name) # Raises error

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

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

84 return dict(self._parameters)

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

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

88 cls = bravais_lattices[self.conventional_cls]

89 return cls(**self._parameters)

91 def tocell(self) -> Cell:

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

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

94 return Cell(cell)

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

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

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

100 # (Just a brute-force implementation)

101 cell = self.tocell()

102 return cell.cellpar()

103

104 @property

105 def special_path(self) -> str:

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

107

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

109 'GXYSGZS1NPY1Z,XP'

110 """

111 return self._variant.special_path

112

113 @property

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

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

116

117 >>> from ase.lattice import BCT

118

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

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

121 """

122 labels = parse_path_string(self._variant.special_point_names)

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

124 return labels[0]

125

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

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

128

129 Ordering is consistent with special_point_names."""

130 if self._variant.special_points is not None:

131 # Fixed dictionary of special points

132 d = self.get_special_points()

133 labels = self.special_point_names

134 assert len(d) == len(labels)

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

136 for i, label in enumerate(labels):

137 points[i] = d[label]

138 return points

139

140 # Special points depend on lattice parameters:

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

142 **self._parameters)

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

144 return np.array(points)

145

146 def get_special_points(self) -> Dict[str, np.ndarray]:

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

148 if self._variant.special_points is not None:

149 return self._variant.special_points

150

151 labels = self.special_point_names

152 points = self.get_special_points_array()

153

154 return dict(zip(labels, points))

155

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

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

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

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

160 npoints=0)

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

162

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

164 density=None) -> BandPath:

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

166

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

168

169 >>> from ase.lattice import BCT

170

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

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

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

174

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

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

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

178 the kpoints may not correspond to your particular cell.

179

180 """

181 if special_points is None:

182 special_points = self.get_special_points()

183

184 if path is None:

185 path = self._variant.special_path

186 elif not isinstance(path, str):

187 from ase.dft.kpoints import resolve_custom_points

188 path, special_points = resolve_custom_points(path,

189 special_points,

190 self._eps)

191

192 cell = self.tocell()

193 bandpath = BandPath(cell=cell, path=path,

194 special_points=special_points)

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

196

197 @abstractmethod

198 def _cell(self, **kwargs):

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

200

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

202

203 def _special_points(self, **kwargs):

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

205

206 Subclasses typically return a nested list.

207

208 Ordering must be same as kpoint labels.

209

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

211 raise NotImplementedError

212

213 def _variant_name(self, **kwargs):

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

215

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

217 raise NotImplementedError

218

219 def __format__(self, spec):

220 tokens = []

221 if not spec:

222 spec = '.6g'

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

224

225 for name in self.parameters:

226 value = self._parameters[name]

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

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

229

230 def __str__(self) -> str:

231 return self.__format__('')

232

233 def __repr__(self) -> str:

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

235

236 def description(self) -> str:

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

238 points = self.get_special_points()

239 labels = self.special_point_names

240

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

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

243 for label in labels])

244

245 string = """\

246{repr}

247 {variant}

248 Special point coordinates:

249{special_points}

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

251 variant=self._variant,

252 special_points=coordstring)

253 return string

254

255 @classmethod

256 def type_description(cls):

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

258 desc = """\

259Lattice name: {name}

260 Long name: {longname}

261 Parameters: {parameters}

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

263

264 chunks = [desc]

265 for name in cls.variant_names:

266 var = cls.variants[name]

267 txt = str(var)

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

269 lines.append('')

270 chunks.extend(lines)

271

272 return '\n'.join(chunks)

273

274

275class Variant:

276 variant_desc = """\

277Variant name: {name}

278 Special point names: {special_point_names}

279 Default path: {special_path}

280"""

281

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

283 special_points=None):

284 self.name = name

285 self.special_point_names = special_point_names

286 self.special_path = special_path

287 if special_points is not None:

288 special_points = dict(special_points)

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

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

291 self.special_points = special_points

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

293 # (easier to transform that way)

294

295 def __str__(self) -> str:

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

297

298

299bravais_names = []

300bravais_lattices = {}

301bravais_classes = {}

302

303

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

305 parameters, variants, ndim=3):

306 """Decorator for Bravais lattice classes.

307

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

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

310

311 def decorate(cls):

312 btype = cls.__name__

313 cls.name = btype

314 cls.longname = longname

315 cls.crystal_family = crystal_family

316 cls.lattice_system = lattice_system

317 cls.pearson_symbol = pearson_symbol

318 cls.parameters = tuple(parameters)

319 cls.variant_names = []

320 cls.variants = {}

321 cls.ndim = ndim

322

323 for [name, special_point_names, special_path,

324 special_points] in variants:

325 cls.variant_names.append(name)

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

327 special_path, special_points)

328

329 # Register in global list and dictionary

330 bravais_names.append(btype)

331 bravais_lattices[btype] = cls

332 bravais_classes[pearson_symbol] = cls

333 return cls

334

335 return decorate

336

337

338# Common mappings of primitive to conventional cells:

339_identity = np.identity(3, int)

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

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

342

343

344class UnconventionalLattice(ValueError):

345 pass

346

347

348class Cubic(BravaisLattice):

349 """Abstract class for cubic lattices."""

350 conventional_cls = 'CUB'

351

352 def __init__(self, a):

353 super().__init__(a=a)

354

355

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

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

358class CUB(Cubic):

359 conventional_cellmap = _identity

360

361 def _cell(self, a):

362 return a * np.eye(3)

363

364

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

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

367class FCC(Cubic):

368 conventional_cellmap = _bcc_map

369

370 def _cell(self, a):

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

372

373

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

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

376class BCC(Cubic):

377 conventional_cellmap = _fcc_map

378

379 def _cell(self, a):

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

381

382

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

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

385 sc_special_points['tetragonal']]])

386class TET(BravaisLattice):

387 conventional_cls = 'TET'

388 conventional_cellmap = _identity

389

390 def __init__(self, a, c):

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

392

393 def _cell(self, a, c):

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

395

396

397# XXX in BCT2 we use S for Sigma.

398# Also in other places I think

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

400 'ac',

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

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

403class BCT(BravaisLattice):

404 conventional_cls = 'TET'

405 conventional_cellmap = _fcc_map

406

407 def __init__(self, a, c):

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

409

410 def _cell(self, a, c):

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

412

413 def _variant_name(self, a, c):

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

415

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

417 a2 = a * a

418 c2 = c * c

419

420 assert variant.name in self.variants

421

422 if variant.name == 'BCT1':

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

424 points = [[0, 0, 0],

425 [-.5, .5, .5],

426 [0., .5, 0.],

427 [.25, .25, .25],

428 [0., 0., .5],

429 [eta, eta, -eta],

430 [-eta, 1 - eta, eta]]

431 else:

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

433 zeta = 0.5 * a2 / c2

434 points = [[0., .0, 0.],

435 [0., .5, 0.],

436 [.25, .25, .25],

437 [-eta, eta, eta],

438 [eta, 1 - eta, -eta],

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

440 [-zeta, zeta, .5],

441 [.5, .5, -zeta],

442 [.5, .5, -.5]]

443 return points

444

445

446def check_orc(a, b, c):

447 if not a < b < c:

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

449 .format(a, b, c))

450

451

452class Orthorhombic(BravaisLattice):

453 """Abstract class for orthorhombic types."""

454

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

456 check_orc(a, b, c)

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

458

459

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

461 'abc',

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

463 sc_special_points['orthorhombic']]])

464class ORC(Orthorhombic):

465 conventional_cls = 'ORC'

466 conventional_cellmap = _identity

467

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

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

470

471

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

473 'oF', 'abc',

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

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

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

477class ORCF(Orthorhombic):

478 conventional_cls = 'ORC'

479 conventional_cellmap = _bcc_map

480

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

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

483

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

485 a2 = a * a

486 b2 = b * b

487 c2 = c * c

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

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

490

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

492 zeta = xminus

493 eta = xplus

494

495 points = [[0, 0, 0],

496 [.5, .5 + zeta, zeta],

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

498 [.5, .5, .5],

499 [1., .5, .5],

500 [0., eta, eta],

501 [1., 1 - eta, 1 - eta],

502 [.5, 0, .5],

503 [.5, .5, 0]]

504 else:

505 assert variant.name == 'ORCF2'

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

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

508 eta = xminus

509

510 points = [[0, 0, 0],

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

512 [.5, .5 + eta, eta],

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

514 [.5 + delta, .5, delta],

515 [.5, .5, .5],

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

517 [phi, .5 + phi, .5],

518 [0., .5, .5],

519 [.5, 0., .5],

520 [.5, .5, 0.]]

521

522 return points

523

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

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

526 if abs(diff) < self._eps:

527 return 'ORCF3'

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

529

530

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

532 'oI', 'abc',

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

534class ORCI(Orthorhombic):

535 conventional_cls = 'ORC'

536 conventional_cellmap = _fcc_map

537

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

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

540

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

542 a2 = a**2

543 b2 = b**2

544 c2 = c**2

545

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

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

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

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

550

551 points = [[0., 0., 0.],

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

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

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

555 [0, .5, 0],

556 [.5, 0, 0],

557 [0., 0., .5],

558 [.25, .25, .25],

559 [-zeta, zeta, zeta],

560 [zeta, 1 - zeta, -zeta],

561 [eta, -eta, eta],

562 [1 - eta, eta, -eta],

563 [.5, .5, -.5]]

564 return points

565

566

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

568 'oC', 'abc',

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

570class ORCC(BravaisLattice):

571 conventional_cls = 'ORC'

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

573

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

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

576 if a >= b:

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

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

579

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

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

582 [0, 0, c]])

583

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

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

586 points = [[0, 0, 0],

587 [zeta, zeta, .5],

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

589 [0, .5, .5],

590 [0, .5, 0],

591 [-.5, .5, .5],

592 [zeta, zeta, 0],

593 [-zeta, 1 - zeta, 0],

594 [-.5, .5, 0],

595 [0, 0, .5]]

596 return points

597

598

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

600 'ac',

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

602 sc_special_points['hexagonal']]])

603class HEX(BravaisLattice):

604 conventional_cls = 'HEX'

605 conventional_cellmap = _identity

606

607 def __init__(self, a, c):

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

609

610 def _cell(self, a, c):

611 x = 0.5 * np.sqrt(3)

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

613 [0., 0., c]])

614

615

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

617 ('a', 'alpha'),

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

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

620class RHL(BravaisLattice):

621 conventional_cls = 'RHL'

622 conventional_cellmap = _identity

623

624 def __init__(self, a, alpha):

625 if alpha >= 120:

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

627 .format(alpha))

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

629

630 def _cell(self, a, alpha):

631 alpha *= np.pi / 180

632 acosa = a * np.cos(alpha)

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

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

635 acosfrac = acosa / acosa2

636 xx = (1 - acosfrac**2)

637 assert xx > 0.0

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

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

640

641 def _variant_name(self, a, alpha):

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

643

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

645 if variant.name == 'RHL1':

646 cosa = np.cos(alpha * _degrees)

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

648 nu = .75 - 0.5 * eta

649 points = [[0, 0, 0],

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

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

652 [.5, .5, 0],

653 [.5, 0, 0],

654 [0, 0, -.5],

655 [eta, nu, nu],

656 [1 - nu, 1 - nu, 1 - eta],

657 [nu, nu, eta - 1],

658 [1 - nu, nu, 0],

659 [nu, 0, -nu],

660 [.5, .5, .5]]

661 else:

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

663 nu = .75 - 0.5 * eta

664 points = [[0, 0, 0],

665 [.5, -.5, 0],

666 [.5, 0, 0],

667 [1 - nu, -nu, 1 - nu],

668 [nu, nu - 1, nu - 1],

669 [eta, eta, eta],

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

671 [.5, -.5, .5]]

672 return points

673

674

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

676 if b > c or alpha >= 90:

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

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

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

680

681

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

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

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

685class MCL(BravaisLattice):

686 conventional_cls = 'MCL'

687 conventional_cellmap = _identity

688

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

690 check_mcl(a, b, c, alpha)

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

692

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

694 alpha *= _degrees

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

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

697

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

699 cosa = np.cos(alpha * _degrees)

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

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

702

703 points = [[0, 0, 0],

704 [.5, .5, 0],

705 [0, .5, .5],

706 [.5, 0, .5],

707 [.5, 0, -.5],

708 [.5, .5, .5],

709 [0, eta, 1 - nu],

710 [0, 1 - eta, nu],

711 [0, eta, -nu],

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

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

714 [.5, eta, -nu],

715 [0, .5, 0],

716 [0, 0, .5],

717 [0, 0, -.5],

718 [.5, 0, 0]]

719 return points

720

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

722 check_mcl(a, b, c, alpha)

723 return 'MCL'

724

725

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

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

728 [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

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

730 ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

731 'GYFLI,I1ZF1,NGM', None],

732 ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

733 'GYFHZIF1,H1Y1XGN,MG', None],

734 ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

735 'GYFHZI,H1Y1XGN,MG', None],

736 ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z',

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

738class MCLC(BravaisLattice):

739 conventional_cls = 'MCL'

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

741

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

743 check_mcl(a, b, c, alpha)

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

745

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

747 alpha *= np.pi / 180

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

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

750

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

752 # from ase.geometry.cell import mclc

753 # okay, this is a bit hacky

754

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

756 # Right now we just repeat the code:

757 check_mcl(a, b, c, alpha)

758

759 a2 = a * a

760 b2 = b * b

761 cosa = np.cos(alpha * _degrees)

762 sina = np.sin(alpha * _degrees)

763 sina2 = sina**2

764

765 cell = self.tocell()

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

767

768 kgamma = lengths_angles[-1]

769

770 eps = self._eps

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

772 # the same precision.

773 if abs(kgamma - 90) < eps:

774 variant = 2

775 elif kgamma > 90:

776 variant = 1

777 elif kgamma < 90:

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

779 if abs(num - 1) < eps:

780 variant = 4

781 elif num < 1:

782 variant = 3

783 else:

784 variant = 5

785 variant = 'MCLC' + str(variant)

786 return variant

787

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

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

790

791 a2 = a * a

792 b2 = b * b

793 # c2 = c * c

794 cosa = np.cos(alpha * _degrees)

795 sina = np.sin(alpha * _degrees)

796 sina2 = sina**2

797

798 if variant == 1 or variant == 2:

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

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

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

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

803

804 points = [[0, 0, 0],

805 [.5, 0, 0],

806 [0, -.5, 0],

807 [1 - zeta, 1 - zeta, 1 - eta],

808 [zeta, zeta, eta],

809 [-zeta, -zeta, 1 - eta],

810 [1 - zeta, -zeta, 1 - eta],

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

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

813 [.5, .5, .5],

814 [.5, 0, .5],

815 [1 - psi, psi - 1, 0],

816 [psi, 1 - psi, 0],

817 [psi - 1, -psi, 0],

818 [.5, .5, 0],

819 [-.5, -.5, 0],

820 [0, 0, .5]]

821 elif variant == 3 or variant == 4:

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

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

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

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

826 phi = 1 + zeta - 2 * mu

827 psi = eta - 2 * delta

828

829 points = [[0, 0, 0],

830 [1 - phi, 1 - phi, 1 - psi],

831 [phi, phi - 1, psi],

832 [1 - phi, -phi, 1 - psi],

833 [zeta, zeta, eta],

834 [1 - zeta, -zeta, 1 - eta],

835 [-zeta, -zeta, 1 - eta],

836 [.5, -.5, .5],

837 [.5, 0, .5],

838 [.5, 0, 0],

839 [0, -.5, 0],

840 [.5, -.5, 0],

841 [mu, mu, delta],

842 [1 - mu, -mu, -delta],

843 [-mu, -mu, -delta],

844 [mu, mu - 1, delta],

845 [0, 0, .5]]

846 elif variant == 5:

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

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

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

850 nu = 2 * mu - zeta

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

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

853 rho = 1 - zeta * a2 / b2

854

855 points = [[0, 0, 0],

856 [nu, nu, omega],

857 [1 - nu, 1 - nu, 1 - omega],

858 [nu, nu - 1, omega],

859 [zeta, zeta, eta],

860 [1 - zeta, -zeta, 1 - eta],

861 [-zeta, -zeta, 1 - eta],

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

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

864 [.5, .5, .5],

865 [.5, 0, .5],

866 [.5, 0, 0],

867 [0, -.5, 0],

868 [.5, -.5, 0],

869 [mu, mu, delta],

870 [1 - mu, -mu, -delta],

871 [-mu, -mu, -delta],

872 [mu, mu - 1, delta],

873 [0, 0, .5]]

874

875 return points

876

877

878tri_angles_explanation = """\

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

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

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

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

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

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

885

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

887

888

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

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

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

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

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

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

895class TRI(BravaisLattice):

896 conventional_cls = 'TRI'

897 conventional_cellmap = _identity

898

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

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

901 gamma=gamma)

902

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

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

905 singamma = np.sin(gamma * _degrees)

906 cosgamma = np.cos(gamma * _degrees)

907 cosbeta = np.cos(beta * _degrees)

908 cosalpha = np.cos(alpha * _degrees)

909 a3x = c * cosbeta

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

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

912 + 2 * cosalpha * cosbeta * cosgamma)

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

914 [a3x, a3y, a3z]])

915

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

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

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

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

920

921 def raise_unconventional():

922 raise UnconventionalLattice(tri_angles_explanation

923 .format(*kangles))

924

925 eps = self._eps

926 if abs(kgamma - 90) < eps:

927 if kalpha > 90 and kbeta > 90:

928 var = '2a'

929 elif kalpha < 90 and kbeta < 90:

930 var = '2b'

931 else:

932 # Is this possible? Maybe due to epsilon

933 raise_unconventional()

934 elif all(kangles > 90):

935 if kgamma > min(kangles):

936 raise_unconventional()

937 var = '1a'

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

939 if kgamma < max(kangles):

940 raise_unconventional()

941 var = '1b'

942 else:

943 raise_unconventional()

944

945 return 'TRI' + var

946

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

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

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

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

951 points = [[0., 0., 0.],

952 [.5, .5, 0],

953 [0, .5, .5],

954 [.5, 0, .5],

955 [.5, .5, .5],

956 [.5, 0, 0],

957 [0, .5, 0],

958 [0, 0, .5]]

959 else:

960 points = [[0, 0, 0],

961 [.5, -.5, 0],

962 [0, 0, .5],

963 [-.5, -.5, .5],

964 [0, -.5, .5],

965 [0, -0.5, 0],

966 [.5, 0, 0],

967 [-.5, 0, .5]]

968

969 return points

970

971

972def get_subset_points(names, points):

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

974 return newpoints

975

976

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

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

979 ndim=2)

980class OBL(BravaisLattice):

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

982 check_rect(a, b)

983 if alpha >= 90:

984 raise UnconventionalLattice(

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

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

987

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

989 cosa = np.cos(alpha * _degrees)

990 sina = np.sin(alpha * _degrees)

991

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

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

994 [0., 0., 0.]])

995

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

997 cosa = np.cos(alpha * _degrees)

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

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

1000

1001 points = [[0, 0, 0],

1002 [0, 0.5, 0],

1003 [eta, 1 - nu, 0],

1004 [.5, .5, 0],

1005 [1 - eta, nu, 0],

1006 [.5, 0, 0]]

1007

1008 return points

1009

1010

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

1012 [['HEX2D', 'GMK', 'GMKG',

1013 get_subset_points('GMK',

1014 sc_special_points['hexagonal'])]],

1015 ndim=2)

1016class HEX2D(BravaisLattice):

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

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

1019

1020 def _cell(self, a):

1021 x = 0.5 * np.sqrt(3)

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

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

1024 [0., 0., 0.]])

1025

1026

1027def check_rect(a, b):

1028 if a >= b:

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

1030

1031

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

1033 [['RECT', 'GXSY', 'GXSYGS',

1034 get_subset_points('GXSY',

1035 sc_special_points['orthorhombic'])]],

1036 ndim=2)

1037class RECT(BravaisLattice):

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

1039 check_rect(a, b)

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

1041

1042 def _cell(self, a, b):

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

1044 [0, b, 0],

1045 [0, 0, 0.]])

1046

1047

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

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

1050class CRECT(BravaisLattice):

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

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

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

1054 # like in ordinary RECT.

1055 #

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

1057 # identical parameters.

1058 if alpha >= 90:

1059 raise UnconventionalLattice(

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

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

1062

1063 def _cell(self, a, alpha):

1064 x = np.cos(alpha * _degrees)

1065 y = np.sin(alpha * _degrees)

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

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

1068 [0, 0, 0.]])

1069

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

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

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

1073 eta = sina2 / sina

1074 cosa = np.cos(alpha * _degrees)

1075 xi = eta * cosa

1076

1077 points = [[0, 0, 0],

1078 [eta, - eta, 0],

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

1080 [0.5, 0.5, 0]]

1081

1082 return points

1083

1084

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

1086 [['SQR', 'GMX', 'MGXM',

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

1088 ndim=2)

1089class SQR(BravaisLattice):

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

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

1092

1093 def _cell(self, a):

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

1095 [0, a, 0],

1096 [0, 0, 0.]])

1097

1098

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

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

1101 ndim=1)

1102class LINE(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, 0.0],

1108 [0.0, 0.0, 0.0],

1109 [0.0, 0.0, 0.0]])

1110

1111

1112def celldiff(cell1, cell2):

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

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

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

1116 v1v2 = cell1.volume * cell2.volume

1117 if v1v2 < 1e-10:

1118 # (Proposed cell may be linearly dependent)

1119 return np.inf

1120

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

1122 x1 = cell1 @ cell1.T

1123 x2 = cell2 @ cell2.T

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

1125 return dev

1126

1127

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

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

1130

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

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

1133

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

1135 lattice, raise RuntimeError."""

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

1137

1138

1139class LatticeMatcher:

1140 def __init__(self, cell, pbc, eps):

1141 self.orig_cell = cell

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

1143 self.cell = cell.uncomplete(pbc)

1144 self.eps = eps

1145 self.niggli_cell, self.niggli_op = self.cell.niggli_reduce(eps=eps)

1146

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

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

1149 # in the table:

1150 self.memory = {}

1151

1152 def match(self, latname, operations):

1153 matches = []

1154

1155 for op_key in operations:

1156 checker_and_op = self.memory.get(op_key)

1157 if checker_and_op is None:

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

1159 # Niggli cell to AFlow form.

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

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

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

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

1164 else:

1165 checker, op_3x3 = checker_and_op

1166

1167 lat, err = checker.query(latname)

1168 if lat is None or err > self.eps:

1169 continue

1170

1171 # This is the full operation encompassing

1172 # both Niggli reduction of user input and mapping the

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

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

1175 matches.append(Match(lat, op))

1176

1177 return matches

1178

1179

1180class Match:

1181 def __init__(self, lat, op):

1182 self.lat = lat

1183 self.op = op

1184

1185 @cached_property

1186 def orthogonality_defect(self):

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

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

1189

1190

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

1192 """Find Bravais lattice representing this cell.

1193

1194 Returns Bravais lattice object representing the cell along with

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

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

1197 from ase.geometry.bravais_type_engine import niggli_op_table

1198

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

1200

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

1202 # the first one we find:

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

1204 # There may be multiple Niggli operations that produce valid

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

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

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

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

1209

1210 if not matches:

1211 continue # Move to next Bravais lattice

1212

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

1214

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

1216 best.op = flip_2d_handedness(best.op)

1217

1218 return best.lat, best.op

1219

1220 raise RuntimeError('Failed to recognize lattice')

1221

1222

1223def flip_2d_handedness(op):

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

1225 # components are necessarily also left-handed which

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

1227 #

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

1229 # z axis and interchanges x/y:

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

1231 return repair_op @ op

1232

1233

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

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

1236# when a match is found.

1237#

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

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

1240lattice_check_orders = {

1241 1: ['LINE'],

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

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

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

1245

1246

1247class NormalizedLatticeMatcher:

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

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

1250

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

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

1253

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

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

1256 by directly through one of the Bravais lattice classes.

1257

1258 Generally for internal use (this module).

1259

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

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

1262 the tolerance eps is not met."""

1263 self.cell = cell

1264 self.eps = eps

1265

1266 self.cellpar = cell.cellpar()

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

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

1269

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

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

1272

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

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

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

1276

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

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

1279 return None

1280

1281 try:

1282 return latcls(*args)

1283 except UnconventionalLattice:

1284 return None

1285

1286 def match(self):

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

1288

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

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

1291 lat, err = self.query(name)

1292 if lat and err < self.eps:

1293 return lat

1294

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

1296 'with lengths and angles {}'

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

1298

1299 def query(self, latname):

1300 """Match cell against named Bravais lattice.

1301

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

1303 meth = getattr(self, latname)

1304

1305 lat = meth()

1306 if lat is None:

1307 return None, None

1308

1309 newcell = lat.tocell()

1310 err = celldiff(self.cell, newcell)

1311 return lat, err

1312

1313 def LINE(self):

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

1315

1316 def SQR(self):

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

1318

1319 def RECT(self):

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

1321

1322 def CRECT(self):

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

1324

1325 def HEX2D(self):

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

1327

1328 def OBL(self):

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

1330

1331 def CUB(self):

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

1333 # it matches, else None.

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

1335

1336 def FCC(self):

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

1338

1339 def BCC(self):

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

1341

1342 def TET(self):

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

1344

1345 def _bct_orci_lengths(self):

1346 # Coordinate-system independent relation for BCT and ORCI

1347 # standard cells:

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

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

1350 # == b² / 2 (ORCI)

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

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

1353 prods = self.prods

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

1355 if any(lengthsqr < 0):

1356 return None

1357 return np.sqrt(lengthsqr)

1358

1359 def BCT(self):

1360 lengths = self._bct_orci_lengths()

1361 if lengths is None:

1362 return None

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

1364

1365 def HEX(self):

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

1367

1368 def RHL(self):

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

1370

1371 def ORC(self):

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

1373

1374 def ORCF(self):

1375 # ORCF standard cell:

1376 # a2 · a3 = a²/4

1377 # a3 · a1 = b²/4

1378 # a1 · a2 = c²/4

1379 prods = self.prods

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

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

1382 return self._check(ORCF, *orcf_abc)

1383

1384 def ORCI(self):

1385 lengths = self._bct_orci_lengths()

1386 if lengths is None:

1387 return None

1388 return self._check(ORCI, *lengths)

1389

1390 def _orcc_ab(self):

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

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

1393 prods = self.prods

1394 orcc_sqr_ab = np.empty(2)

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

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

1397 if all(orcc_sqr_ab > 0):

1398 return np.sqrt(orcc_sqr_ab)

1399

1400 def ORCC(self):

1401 orcc_lengths_ab = self._orcc_ab()

1402 if orcc_lengths_ab is None:

1403 return None

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

1405

1406 def MCL(self):

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

1408

1409 def MCLC(self):

1410 # MCLC is similar to ORCC:

1411 orcc_ab = self._orcc_ab()

1412 if orcc_ab is None:

1413 return None

1414

1415 prods = self.prods

1416 C = self.C

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

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

1419 if -1 < mclc_cosa < 1:

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

1421 if mclc_b > C:

1422 # XXX Temporary fix for certain otherwise

1423 # unrecognizable lattices.

1424 #

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

1426 # something just outside the domain of conventional

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

1428 # propose a nearby conventional lattice instead, which

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

1430 mclc_b = 0.5 * (mclc_b + C)

1431 C = mclc_b

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

1433

1434 def TRI(self):

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

1436

1437

1438def all_variants():

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

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

1441 alpha = 55.0

1442 yield CUB(a)

1443 yield FCC(a)

1444 yield BCC(a)

1445 yield TET(a, c)

1446 bct1 = BCT(2 * a, c)

1447 bct2 = BCT(a, c)

1448 assert bct1.variant == 'BCT1'

1449 assert bct2.variant == 'BCT2'

1450

1451 yield bct1

1452 yield bct2

1453

1454 yield ORC(a, b, c)

1455

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

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

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

1459 orcf3 = ORCF(a0, b, c)

1460 assert orcf1.variant == 'ORCF1'

1461 assert orcf2.variant == 'ORCF2'

1462 assert orcf3.variant == 'ORCF3'

1463 yield orcf1

1464 yield orcf2

1465 yield orcf3

1466

1467 yield ORCI(a, b, c)

1468 yield ORCC(a, b, c)

1469

1470 yield HEX(a, c)

1471

1472 rhl1 = RHL(a, alpha=55.0)

1473 assert rhl1.variant == 'RHL1'

1474 yield rhl1

1475

1476 rhl2 = RHL(a, alpha=105.0)

1477 assert rhl2.variant == 'RHL2'

1478 yield rhl2

1479

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

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

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

1483

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

1485 assert mclc1.variant == 'MCLC1'

1486 yield mclc1

1487 # mclc2 has same special points as mclc1

1488

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

1490 assert mclc3.variant == 'MCLC3'

1491 yield mclc3

1492 # mclc4 has same special points as mclc3

1493

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

1495

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

1497 assert mclc5.variant == 'MCLC5'

1498 yield mclc5

1499

1500 def get_tri(kcellpar):

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

1502 icell = Cell.fromcellpar(kcellpar)

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

1504 return TRI(*cellpar)

1505

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

1507 assert tri1a.variant == 'TRI1a'

1508 yield tri1a

1509

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

1511 assert tri1b.variant == 'TRI1b'

1512 yield tri1b

1513

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

1515 assert tri2a.variant == 'TRI2a'

1516 yield tri2a

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

1518 assert tri2b.variant == 'TRI2b'

1519 yield tri2b

1520

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

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

1523 # of the same lattice.

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

1525 yield RECT(a, b)

1526 yield CRECT(a, alpha=alpha)

1527 yield HEX2D(a)

1528 yield SQR(a)

1529 yield LINE(a)

Coverage for /builds/ase/ase/ase/lattice/init.py: 97.33%

749 statements