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

1from abc import ABC, abstractmethod

2from typing import Dict, List

4import numpy as np

6from ase.cell import Cell

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

8from ase.utils import pbc2pbc

10_degrees = np.pi / 180

13class BravaisLattice(ABC):

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

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

17 five 2D classes.

19 Each class stores basic static information:

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

22 >>> FCC.name

23 'FCC'

24 >>> FCC.longname

25 'face-centred cubic'

26 >>> FCC.pearson_symbol

27 'cF'

28 >>> MCL.parameters

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

31 Each class can be instantiated with the specific lattice parameters

32 that apply to that lattice:

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

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

37 """

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

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

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

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

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

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

44 # 'BCT2': <variant object>}

45 ndim = None

47 def __init__(self, **kwargs):

48 p = {}

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

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

51 p[k] = float(v)

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

53 self._parameters = p

54 self._eps = eps

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

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

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

59 else:

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

61 self._variant = self.variants[name]

63 @property

64 def variant(self) -> str:

65 """Return name of lattice variant.

67 >>> from ase.lattice import BCT

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

70 'BCT2'

71 """

72 return self._variant.name

74 def __getattr__(self, name: str):

75 if name in self._parameters:

76 return self._parameters[name]

77 return self.__getattribute__(name) # Raises error

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

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

81 return dict(self._parameters)

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

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

85 cls = bravais_lattices[self.conventional_cls]

86 return cls(**self._parameters)

88 def tocell(self) -> Cell:

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

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

91 return Cell(cell)

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

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

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

97 # (Just a brute-force implementation)

98 cell = self.tocell()

99 return cell.cellpar()

100

101 @property

102 def special_path(self) -> str:

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

104

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

106 'GXYSGZS1NPY1Z,XP'

107 """

108 return self._variant.special_path

109

110 @property

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

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

113

114 >>> from ase.lattice import BCT

115

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

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

118 """

119 labels = parse_path_string(self._variant.special_point_names)

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

121 return labels[0]

122

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

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

125

126 Ordering is consistent with special_point_names."""

127 if self._variant.special_points is not None:

128 # Fixed dictionary of special points

129 d = self.get_special_points()

130 labels = self.special_point_names

131 assert len(d) == len(labels)

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

133 for i, label in enumerate(labels):

134 points[i] = d[label]

135 return points

136

137 # Special points depend on lattice parameters:

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

139 **self._parameters)

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

141 return np.array(points)

142

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

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

145 if self._variant.special_points is not None:

146 return self._variant.special_points

147

148 labels = self.special_point_names

149 points = self.get_special_points_array()

150

151 return dict(zip(labels, points))

152

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

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

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

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

157 npoints=0)

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

159

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

161 density=None) -> BandPath:

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

163

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

165

166 >>> from ase.lattice import BCT

167

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

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

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

171

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

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

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

175 the kpoints may not correspond to your particular cell.

176

177 """

178 if special_points is None:

179 special_points = self.get_special_points()

180

181 if path is None:

182 path = self._variant.special_path

183 elif not isinstance(path, str):

184 from ase.dft.kpoints import resolve_custom_points

185 path, special_points = resolve_custom_points(path,

186 special_points,

187 self._eps)

188

189 cell = self.tocell()

190 bandpath = BandPath(cell=cell, path=path,

191 special_points=special_points)

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

193

194 @abstractmethod

195 def _cell(self, **kwargs):

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

197

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

199

200 def _special_points(self, **kwargs):

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

202

203 Subclasses typically return a nested list.

204

205 Ordering must be same as kpoint labels.

206

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

208 raise NotImplementedError

209

210 def _variant_name(self, **kwargs):

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

212

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

214 raise NotImplementedError

215

216 def __format__(self, spec):

217 tokens = []

218 if not spec:

219 spec = '.6g'

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

221

222 for name in self.parameters:

223 value = self._parameters[name]

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

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

226

227 def __str__(self) -> str:

228 return self.__format__('')

229

230 def __repr__(self) -> str:

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

232

233 def description(self) -> str:

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

235 points = self.get_special_points()

236 labels = self.special_point_names

237

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

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

240 for label in labels])

241

242 string = """\

243{repr}

244 {variant}

245 Special point coordinates:

246{special_points}

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

248 variant=self._variant,

249 special_points=coordstring)

250 return string

251

252 @classmethod

253 def type_description(cls):

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

255 desc = """\

256Lattice name: {name}

257 Long name: {longname}

258 Parameters: {parameters}

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

260

261 chunks = [desc]

262 for name in cls.variant_names:

263 var = cls.variants[name]

264 txt = str(var)

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

266 lines.append('')

267 chunks.extend(lines)

268

269 return '\n'.join(chunks)

270

271

272class Variant:

273 variant_desc = """\

274Variant name: {name}

275 Special point names: {special_point_names}

276 Default path: {special_path}

277"""

278

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

280 special_points=None):

281 self.name = name

282 self.special_point_names = special_point_names

283 self.special_path = special_path

284 if special_points is not None:

285 special_points = dict(special_points)

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

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

288 self.special_points = special_points

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

290 # (easier to transform that way)

291

292 def __str__(self) -> str:

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

294

295

296bravais_names = []

297bravais_lattices = {}

298bravais_classes = {}

299

300

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

302 parameters, variants, ndim=3):

303 """Decorator for Bravais lattice classes.

304

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

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

307

308 def decorate(cls):

309 btype = cls.__name__

310 cls.name = btype

311 cls.longname = longname

312 cls.crystal_family = crystal_family

313 cls.lattice_system = lattice_system

314 cls.pearson_symbol = pearson_symbol

315 cls.parameters = tuple(parameters)

316 cls.variant_names = []

317 cls.variants = {}

318 cls.ndim = ndim

319

320 for [name, special_point_names, special_path,

321 special_points] in variants:

322 cls.variant_names.append(name)

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

324 special_path, special_points)

325

326 # Register in global list and dictionary

327 bravais_names.append(btype)

328 bravais_lattices[btype] = cls

329 bravais_classes[pearson_symbol] = cls

330 return cls

331

332 return decorate

333

334

335# Common mappings of primitive to conventional cells:

336_identity = np.identity(3, int)

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

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

339

340

341class UnconventionalLattice(ValueError):

342 pass

343

344

345class Cubic(BravaisLattice):

346 """Abstract class for cubic lattices."""

347 conventional_cls = 'CUB'

348

349 def __init__(self, a):

350 super().__init__(a=a)

351

352

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

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

355class CUB(Cubic):

356 conventional_cellmap = _identity

357

358 def _cell(self, a):

359 return a * np.eye(3)

360

361

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

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

364class FCC(Cubic):

365 conventional_cellmap = _bcc_map

366

367 def _cell(self, a):

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

369

370

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

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

373class BCC(Cubic):

374 conventional_cellmap = _fcc_map

375

376 def _cell(self, a):

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

378

379

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

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

382 sc_special_points['tetragonal']]])

383class TET(BravaisLattice):

384 conventional_cls = 'TET'

385 conventional_cellmap = _identity

386

387 def __init__(self, a, c):

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

389

390 def _cell(self, a, c):

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

392

393

394# XXX in BCT2 we use S for Sigma.

395# Also in other places I think

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

397 'ac',

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

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

400class BCT(BravaisLattice):

401 conventional_cls = 'TET'

402 conventional_cellmap = _fcc_map

403

404 def __init__(self, a, c):

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

406

407 def _cell(self, a, c):

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

409

410 def _variant_name(self, a, c):

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

412

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

414 a2 = a * a

415 c2 = c * c

416

417 assert variant.name in self.variants

418

419 if variant.name == 'BCT1':

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

421 points = [[0, 0, 0],

422 [-.5, .5, .5],

423 [0., .5, 0.],

424 [.25, .25, .25],

425 [0., 0., .5],

426 [eta, eta, -eta],

427 [-eta, 1 - eta, eta]]

428 else:

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

430 zeta = 0.5 * a2 / c2

431 points = [[0., .0, 0.],

432 [0., .5, 0.],

433 [.25, .25, .25],

434 [-eta, eta, eta],

435 [eta, 1 - eta, -eta],

436 [0., 0., .5],

437 [-zeta, zeta, .5],

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

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

440 return points

441

442

443def check_orc(a, b, c):

444 if not a < b < c:

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

446 .format(a, b, c))

447

448

449class Orthorhombic(BravaisLattice):

450 """Abstract class for orthorhombic types."""

451

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

453 check_orc(a, b, c)

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

455

456

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

458 'abc',

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

460 sc_special_points['orthorhombic']]])

461class ORC(Orthorhombic):

462 conventional_cls = 'ORC'

463 conventional_cellmap = _identity

464

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

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

467

468

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

470 'oF', 'abc',

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

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

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

474class ORCF(Orthorhombic):

475 conventional_cls = 'ORC'

476 conventional_cellmap = _bcc_map

477

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

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

480

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

482 a2 = a * a

483 b2 = b * b

484 c2 = c * c

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

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

487

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

489 zeta = xminus

490 eta = xplus

491

492 points = [[0, 0, 0],

493 [.5, .5 + zeta, zeta],

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

495 [.5, .5, .5],

496 [1., .5, .5],

497 [0., eta, eta],

498 [1., 1 - eta, 1 - eta],

499 [.5, 0, .5],

500 [.5, .5, 0]]

501 else:

502 assert variant.name == 'ORCF2'

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

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

505 eta = xminus

506

507 points = [[0, 0, 0],

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

509 [.5, .5 + eta, eta],

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

511 [.5 + delta, .5, delta],

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

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

514 [phi, .5 + phi, .5],

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

516 [.5, 0., .5],

517 [.5, .5, 0.]]

518

519 return points

520

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

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

523 if abs(diff) < self._eps:

524 return 'ORCF3'

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

526

527

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

529 'oI', 'abc',

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

531class ORCI(Orthorhombic):

532 conventional_cls = 'ORC'

533 conventional_cellmap = _fcc_map

534

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

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

537

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

539 a2 = a**2

540 b2 = b**2

541 c2 = c**2

542

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

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

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

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

547

548 points = [[0., 0., 0.],

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

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

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

552 [0, .5, 0],

553 [.5, 0, 0],

554 [0., 0., .5],

555 [.25, .25, .25],

556 [-zeta, zeta, zeta],

557 [zeta, 1 - zeta, -zeta],

558 [eta, -eta, eta],

559 [1 - eta, eta, -eta],

560 [.5, .5, -.5]]

561 return points

562

563

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

565 'oC', 'abc',

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

567class ORCC(BravaisLattice):

568 conventional_cls = 'ORC'

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

570

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

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

573 if a >= b:

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

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

576

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

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

579 [0, 0, c]])

580

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

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

583 points = [[0, 0, 0],

584 [zeta, zeta, .5],

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

586 [0, .5, .5],

587 [0, .5, 0],

588 [-.5, .5, .5],

589 [zeta, zeta, 0],

590 [-zeta, 1 - zeta, 0],

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

592 [0, 0, .5]]

593 return points

594

595

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

597 'ac',

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

599 sc_special_points['hexagonal']]])

600class HEX(BravaisLattice):

601 conventional_cls = 'HEX'

602 conventional_cellmap = _identity

603

604 def __init__(self, a, c):

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

606

607 def _cell(self, a, c):

608 x = 0.5 * np.sqrt(3)

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

610 [0., 0., c]])

611

612

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

614 ('a', 'alpha'),

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

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

617class RHL(BravaisLattice):

618 conventional_cls = 'RHL'

619 conventional_cellmap = _identity

620

621 def __init__(self, a, alpha):

622 if alpha >= 120:

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

624 .format(alpha))

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

626

627 def _cell(self, a, alpha):

628 alpha *= np.pi / 180

629 acosa = a * np.cos(alpha)

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

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

632 acosfrac = acosa / acosa2

633 xx = (1 - acosfrac**2)

634 assert xx > 0.0

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

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

637

638 def _variant_name(self, a, alpha):

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

640

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

642 if variant.name == 'RHL1':

643 cosa = np.cos(alpha * _degrees)

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

645 nu = .75 - 0.5 * eta

646 points = [[0, 0, 0],

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

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

649 [.5, .5, 0],

650 [.5, 0, 0],

651 [0, 0, -.5],

652 [eta, nu, nu],

653 [1 - nu, 1 - nu, 1 - eta],

654 [nu, nu, eta - 1],

655 [1 - nu, nu, 0],

656 [nu, 0, -nu],

657 [.5, .5, .5]]

658 else:

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

660 nu = .75 - 0.5 * eta

661 points = [[0, 0, 0],

662 [.5, -.5, 0],

663 [.5, 0, 0],

664 [1 - nu, -nu, 1 - nu],

665 [nu, nu - 1, nu - 1],

666 [eta, eta, eta],

667 [1 - eta, -eta, -eta],

668 [.5, -.5, .5]]

669 return points

670

671

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

673 if b > c or alpha >= 90:

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

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

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

677

678

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

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

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

682class MCL(BravaisLattice):

683 conventional_cls = 'MCL'

684 conventional_cellmap = _identity

685

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

687 check_mcl(a, b, c, alpha)

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

689

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

691 alpha *= _degrees

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

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

694

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

696 cosa = np.cos(alpha * _degrees)

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

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

699

700 points = [[0, 0, 0],

701 [.5, .5, 0],

702 [0, .5, .5],

703 [.5, 0, .5],

704 [.5, 0, -.5],

705 [.5, .5, .5],

706 [0, eta, 1 - nu],

707 [0, 1 - eta, nu],

708 [0, eta, -nu],

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

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

711 [.5, eta, -nu],

712 [0, .5, 0],

713 [0, 0, .5],

714 [0, 0, -.5],

715 [.5, 0, 0]]

716 return points

717

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

719 check_mcl(a, b, c, alpha)

720 return 'MCL'

721

722

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

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

725 [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

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

727 ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

728 'GYFLI,I1ZF1,NGM', None],

729 ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

730 'GYFHZIF1,H1Y1XGN,MG', None],

731 ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

732 'GYFHZI,H1Y1XGN,MG', None],

733 ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z',

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

735class MCLC(BravaisLattice):

736 conventional_cls = 'MCL'

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

738

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

740 check_mcl(a, b, c, alpha)

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

742

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

744 alpha *= np.pi / 180

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

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

747

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

749 # from ase.geometry.cell import mclc

750 # okay, this is a bit hacky

751

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

753 # Right now we just repeat the code:

754 check_mcl(a, b, c, alpha)

755

756 a2 = a * a

757 b2 = b * b

758 cosa = np.cos(alpha * _degrees)

759 sina = np.sin(alpha * _degrees)

760 sina2 = sina**2

761

762 cell = self.tocell()

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

764

765 kgamma = lengths_angles[-1]

766

767 eps = self._eps

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

769 # the same precision.

770 if abs(kgamma - 90) < eps:

771 variant = 2

772 elif kgamma > 90:

773 variant = 1

774 elif kgamma < 90:

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

776 if abs(num - 1) < eps:

777 variant = 4

778 elif num < 1:

779 variant = 3

780 else:

781 variant = 5

782 variant = 'MCLC' + str(variant)

783 return variant

784

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

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

787

788 a2 = a * a

789 b2 = b * b

790 # c2 = c * c

791 cosa = np.cos(alpha * _degrees)

792 sina = np.sin(alpha * _degrees)

793 sina2 = sina**2

794

795 if variant == 1 or variant == 2:

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

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

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

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

800

801 points = [[0, 0, 0],

802 [.5, 0, 0],

803 [0, -.5, 0],

804 [1 - zeta, 1 - zeta, 1 - eta],

805 [zeta, zeta, eta],

806 [-zeta, -zeta, 1 - eta],

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

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

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

810 [.5, .5, .5],

811 [.5, 0, .5],

812 [1 - psi, psi - 1, 0],

813 [psi, 1 - psi, 0],

814 [psi - 1, -psi, 0],

815 [.5, .5, 0],

816 [-.5, -.5, 0],

817 [0, 0, .5]]

818 elif variant == 3 or variant == 4:

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

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

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

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

823 phi = 1 + zeta - 2 * mu

824 psi = eta - 2 * delta

825

826 points = [[0, 0, 0],

827 [1 - phi, 1 - phi, 1 - psi],

828 [phi, phi - 1, psi],

829 [1 - phi, -phi, 1 - psi],

830 [zeta, zeta, eta],

831 [1 - zeta, -zeta, 1 - eta],

832 [-zeta, -zeta, 1 - eta],

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

834 [.5, 0, .5],

835 [.5, 0, 0],

836 [0, -.5, 0],

837 [.5, -.5, 0],

838 [mu, mu, delta],

839 [1 - mu, -mu, -delta],

840 [-mu, -mu, -delta],

841 [mu, mu - 1, delta],

842 [0, 0, .5]]

843 elif variant == 5:

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

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

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

847 nu = 2 * mu - zeta

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

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

850 rho = 1 - zeta * a2 / b2

851

852 points = [[0, 0, 0],

853 [nu, nu, omega],

854 [1 - nu, 1 - nu, 1 - omega],

855 [nu, nu - 1, omega],

856 [zeta, zeta, eta],

857 [1 - zeta, -zeta, 1 - eta],

858 [-zeta, -zeta, 1 - eta],

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

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

861 [.5, .5, .5],

862 [.5, 0, .5],

863 [.5, 0, 0],

864 [0, -.5, 0],

865 [.5, -.5, 0],

866 [mu, mu, delta],

867 [1 - mu, -mu, -delta],

868 [-mu, -mu, -delta],

869 [mu, mu - 1, delta],

870 [0, 0, .5]]

871

872 return points

873

874

875tri_angles_explanation = """\

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

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

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

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

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

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

882

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

884

885

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

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

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

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

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

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

892class TRI(BravaisLattice):

893 conventional_cls = 'TRI'

894 conventional_cellmap = _identity

895

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

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

898 gamma=gamma)

899

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

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

902 singamma = np.sin(gamma * _degrees)

903 cosgamma = np.cos(gamma * _degrees)

904 cosbeta = np.cos(beta * _degrees)

905 cosalpha = np.cos(alpha * _degrees)

906 a3x = c * cosbeta

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

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

909 + 2 * cosalpha * cosbeta * cosgamma)

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

911 [a3x, a3y, a3z]])

912

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

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

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

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

917

918 def raise_unconventional():

919 raise UnconventionalLattice(tri_angles_explanation

920 .format(*kangles))

921

922 eps = self._eps

923 if abs(kgamma - 90) < eps:

924 if kalpha > 90 and kbeta > 90:

925 var = '2a'

926 elif kalpha < 90 and kbeta < 90:

927 var = '2b'

928 else:

929 # Is this possible? Maybe due to epsilon

930 raise_unconventional()

931 elif all(kangles > 90):

932 if kgamma > min(kangles):

933 raise_unconventional()

934 var = '1a'

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

936 if kgamma < max(kangles):

937 raise_unconventional()

938 var = '1b'

939 else:

940 raise_unconventional()

941

942 return 'TRI' + var

943

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

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

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

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

948 points = [[0., 0., 0.],

949 [.5, .5, 0],

950 [0, .5, .5],

951 [.5, 0, .5],

952 [.5, .5, .5],

953 [.5, 0, 0],

954 [0, .5, 0],

955 [0, 0, .5]]

956 else:

957 points = [[0, 0, 0],

958 [.5, -.5, 0],

959 [0, 0, .5],

960 [-.5, -.5, .5],

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

962 [0, -0.5, 0],

963 [.5, 0, 0],

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

965

966 return points

967

968

969def get_subset_points(names, points):

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

971 return newpoints

972

973

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

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

976 ndim=2)

977class OBL(BravaisLattice):

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

979 check_rect(a, b)

980 if alpha >= 90:

981 raise UnconventionalLattice(

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

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

984

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

986 cosa = np.cos(alpha * _degrees)

987 sina = np.sin(alpha * _degrees)

988

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

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

991 [0., 0., 0.]])

992

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

994 cosa = np.cos(alpha * _degrees)

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

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

997

998 points = [[0, 0, 0],

999 [0, 0.5, 0],

1000 [eta, 1 - nu, 0],

1001 [.5, .5, 0],

1002 [1 - eta, nu, 0],

1003 [.5, 0, 0]]

1004

1005 return points

1006

1007

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

1009 [['HEX2D', 'GMK', 'GMKG',

1010 get_subset_points('GMK',

1011 sc_special_points['hexagonal'])]],

1012 ndim=2)

1013class HEX2D(BravaisLattice):

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

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

1016

1017 def _cell(self, a):

1018 x = 0.5 * np.sqrt(3)

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

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

1021 [0., 0., 0.]])

1022

1023

1024def check_rect(a, b):

1025 if a >= b:

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

1027

1028

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

1030 [['RECT', 'GXSY', 'GXSYGS',

1031 get_subset_points('GXSY',

1032 sc_special_points['orthorhombic'])]],

1033 ndim=2)

1034class RECT(BravaisLattice):

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

1036 check_rect(a, b)

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

1038

1039 def _cell(self, a, b):

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

1041 [0, b, 0],

1042 [0, 0, 0.]])

1043

1044

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

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

1047class CRECT(BravaisLattice):

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

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

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

1051 # like in ordinary RECT.

1052 #

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

1054 # identical parameters.

1055 if alpha >= 90:

1056 raise UnconventionalLattice(

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

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

1059

1060 def _cell(self, a, alpha):

1061 x = np.cos(alpha * _degrees)

1062 y = np.sin(alpha * _degrees)

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

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

1065 [0, 0, 0.]])

1066

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

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

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

1070 eta = sina2 / sina

1071 cosa = np.cos(alpha * _degrees)

1072 xi = eta * cosa

1073

1074 points = [[0, 0, 0],

1075 [eta, - eta, 0],

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

1077 [0.5, 0.5, 0]]

1078

1079 return points

1080

1081

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

1083 [['SQR', 'GMX', 'MGXM',

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

1085 ndim=2)

1086class SQR(BravaisLattice):

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

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

1089

1090 def _cell(self, a):

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

1092 [0, a, 0],

1093 [0, 0, 0.]])

1094

1095

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

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

1098 ndim=1)

1099class LINE(BravaisLattice):

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

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

1102

1103 def _cell(self, a):

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

1105 [0.0, 0.0, 0.0],

1106 [0.0, 0.0, 0.0]])

1107

1108

1109def celldiff(cell1, cell2):

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

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

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

1113 v1v2 = cell1.volume * cell2.volume

1114 if v1v2 < 1e-10:

1115 # (Proposed cell may be linearly dependent)

1116 return np.inf

1117

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

1119 x1 = cell1 @ cell1.T

1120 x2 = cell2 @ cell2.T

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

1122 return dev

1123

1124

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

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

1127

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

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

1130

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

1132 lattice, raise RuntimeError."""

1133 return LatticeChecker(cell, eps).match()

1134

1135

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

1137 """Find Bravais lattice representing this cell.

1138

1139 Returns Bravais lattice object representing the cell along with

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

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

1142 from ase.geometry.bravais_type_engine import niggli_op_table

1143

1144 pbc = cell.any(1) & pbc2pbc(pbc)

1145 npbc = sum(pbc)

1146

1147 cell = cell.uncomplete(pbc)

1148 rcell, reduction_op = cell.niggli_reduce(eps=eps)

1149

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

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

1152 # in the table:

1153 memory = {}

1154

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

1156 # the first one we find:

1157 for latname in LatticeChecker.check_orders[npbc]:

1158 # There may be multiple Niggli operations that produce valid

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

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

1161 # just return the first one we find so we must remember then:

1162 matching_lattices = []

1163

1164 for op_key in niggli_op_table[latname]:

1165 checker_and_op = memory.get(op_key)

1166 if checker_and_op is None:

1167 normalization_op = np.array(op_key).reshape(3, 3)

1168 candidate = Cell(np.linalg.inv(normalization_op.T) @ rcell)

1169 checker = LatticeChecker(candidate, eps=eps)

1170 memory[op_key] = (checker, normalization_op)

1171 else:

1172 checker, normalization_op = checker_and_op

1173

1174 lat = checker.query(latname)

1175 if lat is not None:

1176 op = normalization_op @ np.linalg.inv(reduction_op)

1177 matching_lattices.append((lat, op))

1178

1179 # Among any matching lattices, return the one with lowest

1180 # orthogonality defect:

1181 best = None

1182 best_defect = np.inf

1183 for lat, op in matching_lattices:

1184 cell = lat.tocell()

1185 lengths = cell.lengths()[pbc]

1186 generalized_volume = cell.complete().volume

1187 defect = np.prod(lengths) / generalized_volume

1188 if defect < best_defect:

1189 best = lat, op

1190 best_defect = defect

1191

1192 if best is not None:

1193 if npbc == 2:

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

1195 # components are necessarily also left-handed which

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

1197 #

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

1199 # z axis and interchanges x/y:

1200 if op[2, 2] < 0:

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

1202 op = repair_op @ op

1203 best = lat, op

1204

1205 return best

1206

1207 raise RuntimeError('Failed to recognize lattice')

1208

1209

1210class LatticeChecker:

1211 # The check order is slightly different than elsewhere listed order

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

1213 check_orders = {

1214 1: ['LINE'],

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

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

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

1218

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

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

1221

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

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

1224 by directly through one of the Bravais lattice classes.

1225

1226 Generally for internal use (this module).

1227

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

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

1230 the tolerance eps is not met."""

1231 self.cell = cell

1232 self.eps = eps

1233

1234 self.cellpar = cell.cellpar()

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

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

1237

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

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

1240

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

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

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

1244

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

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

1247 return None

1248 try:

1249 lat = latcls(*args)

1250 except UnconventionalLattice:

1251 return None

1252

1253 newcell = lat.tocell()

1254 err = celldiff(self.cell, newcell)

1255 if err < self.eps:

1256 return lat

1257

1258 def match(self):

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

1260

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

1262 for name in self.check_orders[self.cell.rank]:

1263 lat = self.query(name)

1264 if lat:

1265 return lat

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

1267 'with lengths and angles {}'

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

1269

1270 def query(self, latname):

1271 """Match cell against named Bravais lattice.

1272

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

1274 meth = getattr(self, latname)

1275 lat = meth()

1276 return lat

1277

1278 def LINE(self):

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

1280

1281 def SQR(self):

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

1283

1284 def RECT(self):

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

1286

1287 def CRECT(self):

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

1289

1290 def HEX2D(self):

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

1292

1293 def OBL(self):

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

1295

1296 def CUB(self):

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

1298 # it matches, else None.

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

1300

1301 def FCC(self):

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

1303

1304 def BCC(self):

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

1306

1307 def TET(self):

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

1309

1310 def _bct_orci_lengths(self):

1311 # Coordinate-system independent relation for BCT and ORCI

1312 # standard cells:

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

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

1315 # == b² / 2 (ORCI)

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

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

1318 prods = self.prods

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

1320 if any(lengthsqr < 0):

1321 return None

1322 return np.sqrt(lengthsqr)

1323

1324 def BCT(self):

1325 lengths = self._bct_orci_lengths()

1326 if lengths is None:

1327 return None

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

1329

1330 def HEX(self):

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

1332

1333 def RHL(self):

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

1335

1336 def ORC(self):

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

1338

1339 def ORCF(self):

1340 # ORCF standard cell:

1341 # a2 · a3 = a²/4

1342 # a3 · a1 = b²/4

1343 # a1 · a2 = c²/4

1344 prods = self.prods

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

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

1347 return self._check(ORCF, *orcf_abc)

1348

1349 def ORCI(self):

1350 lengths = self._bct_orci_lengths()

1351 if lengths is None:

1352 return None

1353 return self._check(ORCI, *lengths)

1354

1355 def _orcc_ab(self):

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

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

1358 prods = self.prods

1359 orcc_sqr_ab = np.empty(2)

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

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

1362 if all(orcc_sqr_ab > 0):

1363 return np.sqrt(orcc_sqr_ab)

1364

1365 def ORCC(self):

1366 orcc_lengths_ab = self._orcc_ab()

1367 if orcc_lengths_ab is None:

1368 return None

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

1370

1371 def MCL(self):

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

1373

1374 def MCLC(self):

1375 # MCLC is similar to ORCC:

1376 orcc_ab = self._orcc_ab()

1377 if orcc_ab is None:

1378 return None

1379

1380 prods = self.prods

1381 C = self.C

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

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

1384 if -1 < mclc_cosa < 1:

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

1386 if mclc_b > C:

1387 # XXX Temporary fix for certain otherwise

1388 # unrecognizable lattices.

1389 #

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

1391 # something just outside the domain of conventional

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

1393 # propose a nearby conventional lattice instead, which

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

1395 mclc_b = 0.5 * (mclc_b + C)

1396 C = mclc_b

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

1398

1399 def TRI(self):

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

1401

1402

1403def all_variants():

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

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

1406 alpha = 55.0

1407 yield CUB(a)

1408 yield FCC(a)

1409 yield BCC(a)

1410 yield TET(a, c)

1411 bct1 = BCT(2 * a, c)

1412 bct2 = BCT(a, c)

1413 assert bct1.variant == 'BCT1'

1414 assert bct2.variant == 'BCT2'

1415

1416 yield bct1

1417 yield bct2

1418

1419 yield ORC(a, b, c)

1420

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

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

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

1424 orcf3 = ORCF(a0, b, c)

1425 assert orcf1.variant == 'ORCF1'

1426 assert orcf2.variant == 'ORCF2'

1427 assert orcf3.variant == 'ORCF3'

1428 yield orcf1

1429 yield orcf2

1430 yield orcf3

1431

1432 yield ORCI(a, b, c)

1433 yield ORCC(a, b, c)

1434

1435 yield HEX(a, c)

1436

1437 rhl1 = RHL(a, alpha=55.0)

1438 assert rhl1.variant == 'RHL1'

1439 yield rhl1

1440

1441 rhl2 = RHL(a, alpha=105.0)

1442 assert rhl2.variant == 'RHL2'

1443 yield rhl2

1444

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

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

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

1448

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

1450 assert mclc1.variant == 'MCLC1'

1451 yield mclc1

1452 # mclc2 has same special points as mclc1

1453

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

1455 assert mclc3.variant == 'MCLC3'

1456 yield mclc3

1457 # mclc4 has same special points as mclc3

1458

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

1460

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

1462 assert mclc5.variant == 'MCLC5'

1463 yield mclc5

1464

1465 def get_tri(kcellpar):

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

1467 icell = Cell.fromcellpar(kcellpar)

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

1469 return TRI(*cellpar)

1470

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

1472 assert tri1a.variant == 'TRI1a'

1473 yield tri1a

1474

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

1476 assert tri1b.variant == 'TRI1b'

1477 yield tri1b

1478

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

1480 assert tri2a.variant == 'TRI2a'

1481 yield tri2a

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

1483 assert tri2b.variant == 'TRI2b'

1484 yield tri2b

1485

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

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

1488 # of the same lattice.

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

1490 yield RECT(a, b)

1491 yield CRECT(a, alpha=alpha)

1492 yield HEX2D(a)

1493 yield SQR(a)

1494 yield LINE(a)

Coverage for /builds/hweiske/ase/ase/lattice/init.py: 97.30%

740 statements