Coverage for /builds/kinetik161/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 not (b <= c and 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 = {}

971 for name in names:

972 newpoints[name] = points[name]

973

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 LatticeChecker(cell, eps).match()

1137

1138

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

1140 """Find Bravais lattice representing this cell.

1141

1142 Returns Bravais lattice object representing the cell along with

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

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

1145 from ase.geometry.bravais_type_engine import niggli_op_table

1146

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

1148 npbc = sum(pbc)

1149

1150 cell = cell.uncomplete(pbc)

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

1152

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

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

1155 # in the table:

1156 memory = {}

1157

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

1159 # the first one we find:

1160 for latname in LatticeChecker.check_orders[npbc]:

1161 # There may be multiple Niggli operations that produce valid

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

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

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

1165 matching_lattices = []

1166

1167 for op_key in niggli_op_table[latname]:

1168 checker_and_op = memory.get(op_key)

1169 if checker_and_op is None:

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

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

1172 checker = LatticeChecker(candidate, eps=eps)

1173 memory[op_key] = (checker, normalization_op)

1174 else:

1175 checker, normalization_op = checker_and_op

1176

1177 lat = checker.query(latname)

1178 if lat is not None:

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

1180 matching_lattices.append((lat, op))

1181

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

1183 # orthogonality defect:

1184 best = None

1185 best_defect = np.inf

1186 for lat, op in matching_lattices:

1187 cell = lat.tocell()

1188 lengths = cell.lengths()[pbc]

1189 generalized_volume = cell.complete().volume

1190 defect = np.prod(lengths) / generalized_volume

1191 if defect < best_defect:

1192 best = lat, op

1193 best_defect = defect

1194

1195 if best is not None:

1196 if npbc == 2:

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

1198 # components are necessarily also left-handed which

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

1200 #

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

1202 # z axis and interchanges x/y:

1203 if op[2, 2] < 0:

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

1205 op = repair_op @ op

1206 best = lat, op

1207

1208 return best

1209

1210 raise RuntimeError('Failed to recognize lattice')

1211

1212

1213class LatticeChecker:

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

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

1216 check_orders = {

1217 1: ['LINE'],

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

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

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

1221

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

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

1224

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

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

1227 by directly through one of the Bravais lattice classes.

1228

1229 Generally for internal use (this module).

1230

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

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

1233 the tolerance eps is not met."""

1234 self.cell = cell

1235 self.eps = eps

1236

1237 self.cellpar = cell.cellpar()

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

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

1240

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

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

1243

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

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

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

1247

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

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

1250 return None

1251 try:

1252 lat = latcls(*args)

1253 except UnconventionalLattice:

1254 return None

1255

1256 newcell = lat.tocell()

1257 err = celldiff(self.cell, newcell)

1258 if err < self.eps:

1259 return lat

1260

1261 def match(self):

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

1263

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

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

1266 lat = self.query(name)

1267 if lat:

1268 return lat

1269 else:

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

1271 'with lengths and angles {}'

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

1273

1274 def query(self, latname):

1275 """Match cell against named Bravais lattice.

1276

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

1278 meth = getattr(self, latname)

1279 lat = meth()

1280 return lat

1281

1282 def LINE(self):

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

1284

1285 def SQR(self):

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

1287

1288 def RECT(self):

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

1290

1291 def CRECT(self):

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

1293

1294 def HEX2D(self):

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

1296

1297 def OBL(self):

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

1299

1300 def CUB(self):

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

1302 # it matches, else None.

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

1304

1305 def FCC(self):

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

1307

1308 def BCC(self):

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

1310

1311 def TET(self):

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

1313

1314 def _bct_orci_lengths(self):

1315 # Coordinate-system independent relation for BCT and ORCI

1316 # standard cells:

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

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

1319 # == b² / 2 (ORCI)

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

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

1322 prods = self.prods

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

1324 if any(lengthsqr < 0):

1325 return None

1326 return np.sqrt(lengthsqr)

1327

1328 def BCT(self):

1329 lengths = self._bct_orci_lengths()

1330 if lengths is None:

1331 return None

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

1333

1334 def HEX(self):

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

1336

1337 def RHL(self):

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

1339

1340 def ORC(self):

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

1342

1343 def ORCF(self):

1344 # ORCF standard cell:

1345 # a2 · a3 = a²/4

1346 # a3 · a1 = b²/4

1347 # a1 · a2 = c²/4

1348 prods = self.prods

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

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

1351 return self._check(ORCF, *orcf_abc)

1352

1353 def ORCI(self):

1354 lengths = self._bct_orci_lengths()

1355 if lengths is None:

1356 return None

1357 return self._check(ORCI, *lengths)

1358

1359 def _orcc_ab(self):

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

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

1362 prods = self.prods

1363 orcc_sqr_ab = np.empty(2)

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

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

1366 if all(orcc_sqr_ab > 0):

1367 return np.sqrt(orcc_sqr_ab)

1368

1369 def ORCC(self):

1370 orcc_lengths_ab = self._orcc_ab()

1371 if orcc_lengths_ab is None:

1372 return None

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

1374

1375 def MCL(self):

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

1377

1378 def MCLC(self):

1379 # MCLC is similar to ORCC:

1380 orcc_ab = self._orcc_ab()

1381 if orcc_ab is None:

1382 return None

1383

1384 prods = self.prods

1385 C = self.C

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

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

1388 if -1 < mclc_cosa < 1:

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

1390 if mclc_b > C:

1391 # XXX Temporary fix for certain otherwise

1392 # unrecognizable lattices.

1393 #

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

1395 # something just outside the domain of conventional

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

1397 # propose a nearby conventional lattice instead, which

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

1399 mclc_b = 0.5 * (mclc_b + C)

1400 C = mclc_b

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

1402

1403 def TRI(self):

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

1405

1406

1407def all_variants():

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

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

1410 alpha = 55.0

1411 yield CUB(a)

1412 yield FCC(a)

1413 yield BCC(a)

1414 yield TET(a, c)

1415 bct1 = BCT(2 * a, c)

1416 bct2 = BCT(a, c)

1417 assert bct1.variant == 'BCT1'

1418 assert bct2.variant == 'BCT2'

1419

1420 yield bct1

1421 yield bct2

1422

1423 yield ORC(a, b, c)

1424

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

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

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

1428 orcf3 = ORCF(a0, b, c)

1429 assert orcf1.variant == 'ORCF1'

1430 assert orcf2.variant == 'ORCF2'

1431 assert orcf3.variant == 'ORCF3'

1432 yield orcf1

1433 yield orcf2

1434 yield orcf3

1435

1436 yield ORCI(a, b, c)

1437 yield ORCC(a, b, c)

1438

1439 yield HEX(a, c)

1440

1441 rhl1 = RHL(a, alpha=55.0)

1442 assert rhl1.variant == 'RHL1'

1443 yield rhl1

1444

1445 rhl2 = RHL(a, alpha=105.0)

1446 assert rhl2.variant == 'RHL2'

1447 yield rhl2

1448

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

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

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

1452

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

1454 assert mclc1.variant == 'MCLC1'

1455 yield mclc1

1456 # mclc2 has same special points as mclc1

1457

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

1459 assert mclc3.variant == 'MCLC3'

1460 yield mclc3

1461 # mclc4 has same special points as mclc3

1462

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

1464

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

1466 assert mclc5.variant == 'MCLC5'

1467 yield mclc5

1468

1469 def get_tri(kcellpar):

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

1471 icell = Cell.fromcellpar(kcellpar)

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

1473 return TRI(*cellpar)

1474

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

1476 assert tri1a.variant == 'TRI1a'

1477 yield tri1a

1478

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

1480 assert tri1b.variant == 'TRI1b'

1481 yield tri1b

1482

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

1484 assert tri2a.variant == 'TRI2a'

1485 yield tri2a

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

1487 assert tri2b.variant == 'TRI2b'

1488 yield tri2b

1489

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

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

1492 # of the same lattice.

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

1494 yield RECT(a, b)

1495 yield CRECT(a, alpha=alpha)

1496 yield HEX2D(a)

1497 yield SQR(a)

1498 yield LINE(a)

Coverage for /builds/kinetik161/ase/ase/lattice/init.py: 97.04%

742 statements