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

1# fmt: off

3from abc import ABC, abstractmethod

4from typing import Dict, List

6import numpy as np

8from ase.cell import Cell

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

10from ase.utils import pbc2pbc

12_degrees = np.pi / 180

15class BravaisLattice(ABC):

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

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

19 five 2D classes.

21 Each class stores basic static information:

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

24 >>> FCC.name

25 'FCC'

26 >>> FCC.longname

27 'face-centred cubic'

28 >>> FCC.pearson_symbol

29 'cF'

30 >>> MCL.parameters

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

33 Each class can be instantiated with the specific lattice parameters

34 that apply to that lattice:

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

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

39 """

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

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

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

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

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

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

46 # 'BCT2': <variant object>}

47 ndim = None

49 def __init__(self, **kwargs):

50 p = {}

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

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

53 p[k] = float(v)

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

55 self._parameters = p

56 self._eps = eps

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

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

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

61 else:

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

63 self._variant = self.variants[name]

65 @property

66 def variant(self) -> str:

67 """Return name of lattice variant.

69 >>> from ase.lattice import BCT

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

72 'BCT2'

73 """

74 return self._variant.name

76 def __getattr__(self, name: str):

77 if name in self._parameters:

78 return self._parameters[name]

79 return self.__getattribute__(name) # Raises error

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

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

83 return dict(self._parameters)

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

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

87 cls = bravais_lattices[self.conventional_cls]

88 return cls(**self._parameters)

90 def tocell(self) -> Cell:

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

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

93 return Cell(cell)

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

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

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

99 # (Just a brute-force implementation)

100 cell = self.tocell()

101 return cell.cellpar()

102

103 @property

104 def special_path(self) -> str:

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

106

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

108 'GXYSGZS1NPY1Z,XP'

109 """

110 return self._variant.special_path

111

112 @property

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

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

115

116 >>> from ase.lattice import BCT

117

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

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

120 """

121 labels = parse_path_string(self._variant.special_point_names)

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

123 return labels[0]

124

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

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

127

128 Ordering is consistent with special_point_names."""

129 if self._variant.special_points is not None:

130 # Fixed dictionary of special points

131 d = self.get_special_points()

132 labels = self.special_point_names

133 assert len(d) == len(labels)

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

135 for i, label in enumerate(labels):

136 points[i] = d[label]

137 return points

138

139 # Special points depend on lattice parameters:

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

141 **self._parameters)

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

143 return np.array(points)

144

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

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

147 if self._variant.special_points is not None:

148 return self._variant.special_points

149

150 labels = self.special_point_names

151 points = self.get_special_points_array()

152

153 return dict(zip(labels, points))

154

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

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

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

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

159 npoints=0)

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

161

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

163 density=None) -> BandPath:

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

165

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

167

168 >>> from ase.lattice import BCT

169

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

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

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

173

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

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

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

177 the kpoints may not correspond to your particular cell.

178

179 """

180 if special_points is None:

181 special_points = self.get_special_points()

182

183 if path is None:

184 path = self._variant.special_path

185 elif not isinstance(path, str):

186 from ase.dft.kpoints import resolve_custom_points

187 path, special_points = resolve_custom_points(path,

188 special_points,

189 self._eps)

190

191 cell = self.tocell()

192 bandpath = BandPath(cell=cell, path=path,

193 special_points=special_points)

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

195

196 @abstractmethod

197 def _cell(self, **kwargs):

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

199

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

201

202 def _special_points(self, **kwargs):

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

204

205 Subclasses typically return a nested list.

206

207 Ordering must be same as kpoint labels.

208

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

210 raise NotImplementedError

211

212 def _variant_name(self, **kwargs):

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

214

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

216 raise NotImplementedError

217

218 def __format__(self, spec):

219 tokens = []

220 if not spec:

221 spec = '.6g'

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

223

224 for name in self.parameters:

225 value = self._parameters[name]

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

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

228

229 def __str__(self) -> str:

230 return self.__format__('')

231

232 def __repr__(self) -> str:

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

234

235 def description(self) -> str:

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

237 points = self.get_special_points()

238 labels = self.special_point_names

239

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

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

242 for label in labels])

243

244 string = """\

245{repr}

246 {variant}

247 Special point coordinates:

248{special_points}

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

250 variant=self._variant,

251 special_points=coordstring)

252 return string

253

254 @classmethod

255 def type_description(cls):

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

257 desc = """\

258Lattice name: {name}

259 Long name: {longname}

260 Parameters: {parameters}

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

262

263 chunks = [desc]

264 for name in cls.variant_names:

265 var = cls.variants[name]

266 txt = str(var)

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

268 lines.append('')

269 chunks.extend(lines)

270

271 return '\n'.join(chunks)

272

273

274class Variant:

275 variant_desc = """\

276Variant name: {name}

277 Special point names: {special_point_names}

278 Default path: {special_path}

279"""

280

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

282 special_points=None):

283 self.name = name

284 self.special_point_names = special_point_names

285 self.special_path = special_path

286 if special_points is not None:

287 special_points = dict(special_points)

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

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

290 self.special_points = special_points

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

292 # (easier to transform that way)

293

294 def __str__(self) -> str:

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

296

297

298bravais_names = []

299bravais_lattices = {}

300bravais_classes = {}

301

302

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

304 parameters, variants, ndim=3):

305 """Decorator for Bravais lattice classes.

306

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

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

309

310 def decorate(cls):

311 btype = cls.__name__

312 cls.name = btype

313 cls.longname = longname

314 cls.crystal_family = crystal_family

315 cls.lattice_system = lattice_system

316 cls.pearson_symbol = pearson_symbol

317 cls.parameters = tuple(parameters)

318 cls.variant_names = []

319 cls.variants = {}

320 cls.ndim = ndim

321

322 for [name, special_point_names, special_path,

323 special_points] in variants:

324 cls.variant_names.append(name)

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

326 special_path, special_points)

327

328 # Register in global list and dictionary

329 bravais_names.append(btype)

330 bravais_lattices[btype] = cls

331 bravais_classes[pearson_symbol] = cls

332 return cls

333

334 return decorate

335

336

337# Common mappings of primitive to conventional cells:

338_identity = np.identity(3, int)

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

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

341

342

343class UnconventionalLattice(ValueError):

344 pass

345

346

347class Cubic(BravaisLattice):

348 """Abstract class for cubic lattices."""

349 conventional_cls = 'CUB'

350

351 def __init__(self, a):

352 super().__init__(a=a)

353

354

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

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

357class CUB(Cubic):

358 conventional_cellmap = _identity

359

360 def _cell(self, a):

361 return a * np.eye(3)

362

363

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

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

366class FCC(Cubic):

367 conventional_cellmap = _bcc_map

368

369 def _cell(self, a):

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

371

372

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

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

375class BCC(Cubic):

376 conventional_cellmap = _fcc_map

377

378 def _cell(self, a):

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

380

381

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

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

384 sc_special_points['tetragonal']]])

385class TET(BravaisLattice):

386 conventional_cls = 'TET'

387 conventional_cellmap = _identity

388

389 def __init__(self, a, c):

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

391

392 def _cell(self, a, c):

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

394

395

396# XXX in BCT2 we use S for Sigma.

397# Also in other places I think

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

399 'ac',

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

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

402class BCT(BravaisLattice):

403 conventional_cls = 'TET'

404 conventional_cellmap = _fcc_map

405

406 def __init__(self, a, c):

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

408

409 def _cell(self, a, c):

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

411

412 def _variant_name(self, a, c):

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

414

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

416 a2 = a * a

417 c2 = c * c

418

419 assert variant.name in self.variants

420

421 if variant.name == 'BCT1':

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

423 points = [[0, 0, 0],

424 [-.5, .5, .5],

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

426 [.25, .25, .25],

427 [0., 0., .5],

428 [eta, eta, -eta],

429 [-eta, 1 - eta, eta]]

430 else:

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

432 zeta = 0.5 * a2 / c2

433 points = [[0., .0, 0.],

434 [0., .5, 0.],

435 [.25, .25, .25],

436 [-eta, eta, eta],

437 [eta, 1 - eta, -eta],

438 [0., 0., .5],

439 [-zeta, zeta, .5],

440 [.5, .5, -zeta],

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

442 return points

443

444

445def check_orc(a, b, c):

446 if not a < b < c:

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

448 .format(a, b, c))

449

450

451class Orthorhombic(BravaisLattice):

452 """Abstract class for orthorhombic types."""

453

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

455 check_orc(a, b, c)

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

457

458

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

460 'abc',

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

462 sc_special_points['orthorhombic']]])

463class ORC(Orthorhombic):

464 conventional_cls = 'ORC'

465 conventional_cellmap = _identity

466

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

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

469

470

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

472 'oF', 'abc',

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

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

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

476class ORCF(Orthorhombic):

477 conventional_cls = 'ORC'

478 conventional_cellmap = _bcc_map

479

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

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

482

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

484 a2 = a * a

485 b2 = b * b

486 c2 = c * c

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

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

489

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

491 zeta = xminus

492 eta = xplus

493

494 points = [[0, 0, 0],

495 [.5, .5 + zeta, zeta],

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

497 [.5, .5, .5],

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

499 [0., eta, eta],

500 [1., 1 - eta, 1 - eta],

501 [.5, 0, .5],

502 [.5, .5, 0]]

503 else:

504 assert variant.name == 'ORCF2'

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

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

507 eta = xminus

508

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

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

511 [.5, .5 + eta, eta],

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

513 [.5 + delta, .5, delta],

514 [.5, .5, .5],

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

516 [phi, .5 + phi, .5],

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

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

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

520

521 return points

522

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

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

525 if abs(diff) < self._eps:

526 return 'ORCF3'

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

528

529

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

531 'oI', 'abc',

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

533class ORCI(Orthorhombic):

534 conventional_cls = 'ORC'

535 conventional_cellmap = _fcc_map

536

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

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

539

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

541 a2 = a**2

542 b2 = b**2

543 c2 = c**2

544

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

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

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

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

549

550 points = [[0., 0., 0.],

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

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

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

554 [0, .5, 0],

555 [.5, 0, 0],

556 [0., 0., .5],

557 [.25, .25, .25],

558 [-zeta, zeta, zeta],

559 [zeta, 1 - zeta, -zeta],

560 [eta, -eta, eta],

561 [1 - eta, eta, -eta],

562 [.5, .5, -.5]]

563 return points

564

565

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

567 'oC', 'abc',

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

569class ORCC(BravaisLattice):

570 conventional_cls = 'ORC'

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

572

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

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

575 if a >= b:

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

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

578

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

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

581 [0, 0, c]])

582

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

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

585 points = [[0, 0, 0],

586 [zeta, zeta, .5],

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

588 [0, .5, .5],

589 [0, .5, 0],

590 [-.5, .5, .5],

591 [zeta, zeta, 0],

592 [-zeta, 1 - zeta, 0],

593 [-.5, .5, 0],

594 [0, 0, .5]]

595 return points

596

597

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

599 'ac',

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

601 sc_special_points['hexagonal']]])

602class HEX(BravaisLattice):

603 conventional_cls = 'HEX'

604 conventional_cellmap = _identity

605

606 def __init__(self, a, c):

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

608

609 def _cell(self, a, c):

610 x = 0.5 * np.sqrt(3)

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

612 [0., 0., c]])

613

614

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

616 ('a', 'alpha'),

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

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

619class RHL(BravaisLattice):

620 conventional_cls = 'RHL'

621 conventional_cellmap = _identity

622

623 def __init__(self, a, alpha):

624 if alpha >= 120:

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

626 .format(alpha))

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

628

629 def _cell(self, a, alpha):

630 alpha *= np.pi / 180

631 acosa = a * np.cos(alpha)

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

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

634 acosfrac = acosa / acosa2

635 xx = (1 - acosfrac**2)

636 assert xx > 0.0

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

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

639

640 def _variant_name(self, a, alpha):

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

642

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

644 if variant.name == 'RHL1':

645 cosa = np.cos(alpha * _degrees)

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

647 nu = .75 - 0.5 * eta

648 points = [[0, 0, 0],

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

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

651 [.5, .5, 0],

652 [.5, 0, 0],

653 [0, 0, -.5],

654 [eta, nu, nu],

655 [1 - nu, 1 - nu, 1 - eta],

656 [nu, nu, eta - 1],

657 [1 - nu, nu, 0],

658 [nu, 0, -nu],

659 [.5, .5, .5]]

660 else:

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

662 nu = .75 - 0.5 * eta

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

664 [.5, -.5, 0],

665 [.5, 0, 0],

666 [1 - nu, -nu, 1 - nu],

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

668 [eta, eta, eta],

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

670 [.5, -.5, .5]]

671 return points

672

673

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

675 if b > c or alpha >= 90:

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

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

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

679

680

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

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

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

684class MCL(BravaisLattice):

685 conventional_cls = 'MCL'

686 conventional_cellmap = _identity

687

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

689 check_mcl(a, b, c, alpha)

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

691

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

693 alpha *= _degrees

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

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

696

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

698 cosa = np.cos(alpha * _degrees)

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

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

701

702 points = [[0, 0, 0],

703 [.5, .5, 0],

704 [0, .5, .5],

705 [.5, 0, .5],

706 [.5, 0, -.5],

707 [.5, .5, .5],

708 [0, eta, 1 - nu],

709 [0, 1 - eta, nu],

710 [0, eta, -nu],

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

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

713 [.5, eta, -nu],

714 [0, .5, 0],

715 [0, 0, .5],

716 [0, 0, -.5],

717 [.5, 0, 0]]

718 return points

719

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

721 check_mcl(a, b, c, alpha)

722 return 'MCL'

723

724

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

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

727 [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

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

729 ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z',

730 'GYFLI,I1ZF1,NGM', None],

731 ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

732 'GYFHZIF1,H1Y1XGN,MG', None],

733 ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',

734 'GYFHZI,H1Y1XGN,MG', None],

735 ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z',

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

737class MCLC(BravaisLattice):

738 conventional_cls = 'MCL'

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

740

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

742 check_mcl(a, b, c, alpha)

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

744

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

746 alpha *= np.pi / 180

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

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

749

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

751 # from ase.geometry.cell import mclc

752 # okay, this is a bit hacky

753

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

755 # Right now we just repeat the code:

756 check_mcl(a, b, c, alpha)

757

758 a2 = a * a

759 b2 = b * b

760 cosa = np.cos(alpha * _degrees)

761 sina = np.sin(alpha * _degrees)

762 sina2 = sina**2

763

764 cell = self.tocell()

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

766

767 kgamma = lengths_angles[-1]

768

769 eps = self._eps

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

771 # the same precision.

772 if abs(kgamma - 90) < eps:

773 variant = 2

774 elif kgamma > 90:

775 variant = 1

776 elif kgamma < 90:

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

778 if abs(num - 1) < eps:

779 variant = 4

780 elif num < 1:

781 variant = 3

782 else:

783 variant = 5

784 variant = 'MCLC' + str(variant)

785 return variant

786

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

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

789

790 a2 = a * a

791 b2 = b * b

792 # c2 = c * c

793 cosa = np.cos(alpha * _degrees)

794 sina = np.sin(alpha * _degrees)

795 sina2 = sina**2

796

797 if variant == 1 or variant == 2:

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

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

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

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

802

803 points = [[0, 0, 0],

804 [.5, 0, 0],

805 [0, -.5, 0],

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

807 [zeta, zeta, eta],

808 [-zeta, -zeta, 1 - eta],

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

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

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

812 [.5, .5, .5],

813 [.5, 0, .5],

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

815 [psi, 1 - psi, 0],

816 [psi - 1, -psi, 0],

817 [.5, .5, 0],

818 [-.5, -.5, 0],

819 [0, 0, .5]]

820 elif variant == 3 or variant == 4:

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

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

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

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

825 phi = 1 + zeta - 2 * mu

826 psi = eta - 2 * delta

827

828 points = [[0, 0, 0],

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

830 [phi, phi - 1, psi],

831 [1 - phi, -phi, 1 - psi],

832 [zeta, zeta, eta],

833 [1 - zeta, -zeta, 1 - eta],

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

835 [.5, -.5, .5],

836 [.5, 0, .5],

837 [.5, 0, 0],

838 [0, -.5, 0],

839 [.5, -.5, 0],

840 [mu, mu, delta],

841 [1 - mu, -mu, -delta],

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

843 [mu, mu - 1, delta],

844 [0, 0, .5]]

845 elif variant == 5:

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

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

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

849 nu = 2 * mu - zeta

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

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

852 rho = 1 - zeta * a2 / b2

853

854 points = [[0, 0, 0],

855 [nu, nu, omega],

856 [1 - nu, 1 - nu, 1 - omega],

857 [nu, nu - 1, omega],

858 [zeta, zeta, eta],

859 [1 - zeta, -zeta, 1 - eta],

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

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

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

863 [.5, .5, .5],

864 [.5, 0, .5],

865 [.5, 0, 0],

866 [0, -.5, 0],

867 [.5, -.5, 0],

868 [mu, mu, delta],

869 [1 - mu, -mu, -delta],

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

871 [mu, mu - 1, delta],

872 [0, 0, .5]]

873

874 return points

875

876

877tri_angles_explanation = """\

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

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

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

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

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

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

884

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

886

887

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

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

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

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

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

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

894class TRI(BravaisLattice):

895 conventional_cls = 'TRI'

896 conventional_cellmap = _identity

897

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

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

900 gamma=gamma)

901

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

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

904 singamma = np.sin(gamma * _degrees)

905 cosgamma = np.cos(gamma * _degrees)

906 cosbeta = np.cos(beta * _degrees)

907 cosalpha = np.cos(alpha * _degrees)

908 a3x = c * cosbeta

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

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

911 + 2 * cosalpha * cosbeta * cosgamma)

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

913 [a3x, a3y, a3z]])

914

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

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

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

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

919

920 def raise_unconventional():

921 raise UnconventionalLattice(tri_angles_explanation

922 .format(*kangles))

923

924 eps = self._eps

925 if abs(kgamma - 90) < eps:

926 if kalpha > 90 and kbeta > 90:

927 var = '2a'

928 elif kalpha < 90 and kbeta < 90:

929 var = '2b'

930 else:

931 # Is this possible? Maybe due to epsilon

932 raise_unconventional()

933 elif all(kangles > 90):

934 if kgamma > min(kangles):

935 raise_unconventional()

936 var = '1a'

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

938 if kgamma < max(kangles):

939 raise_unconventional()

940 var = '1b'

941 else:

942 raise_unconventional()

943

944 return 'TRI' + var

945

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

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

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

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

950 points = [[0., 0., 0.],

951 [.5, .5, 0],

952 [0, .5, .5],

953 [.5, 0, .5],

954 [.5, .5, .5],

955 [.5, 0, 0],

956 [0, .5, 0],

957 [0, 0, .5]]

958 else:

959 points = [[0, 0, 0],

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

961 [0, 0, .5],

962 [-.5, -.5, .5],

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

964 [0, -0.5, 0],

965 [.5, 0, 0],

966 [-.5, 0, .5]]

967

968 return points

969

970

971def get_subset_points(names, points):

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

973 return newpoints

974

975

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

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

978 ndim=2)

979class OBL(BravaisLattice):

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

981 check_rect(a, b)

982 if alpha >= 90:

983 raise UnconventionalLattice(

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

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

986

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

988 cosa = np.cos(alpha * _degrees)

989 sina = np.sin(alpha * _degrees)

990

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

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

993 [0., 0., 0.]])

994

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

996 cosa = np.cos(alpha * _degrees)

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

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

999

1000 points = [[0, 0, 0],

1001 [0, 0.5, 0],

1002 [eta, 1 - nu, 0],

1003 [.5, .5, 0],

1004 [1 - eta, nu, 0],

1005 [.5, 0, 0]]

1006

1007 return points

1008

1009

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

1011 [['HEX2D', 'GMK', 'GMKG',

1012 get_subset_points('GMK',

1013 sc_special_points['hexagonal'])]],

1014 ndim=2)

1015class HEX2D(BravaisLattice):

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

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

1018

1019 def _cell(self, a):

1020 x = 0.5 * np.sqrt(3)

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

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

1023 [0., 0., 0.]])

1024

1025

1026def check_rect(a, b):

1027 if a >= b:

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

1029

1030

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

1032 [['RECT', 'GXSY', 'GXSYGS',

1033 get_subset_points('GXSY',

1034 sc_special_points['orthorhombic'])]],

1035 ndim=2)

1036class RECT(BravaisLattice):

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

1038 check_rect(a, b)

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

1040

1041 def _cell(self, a, b):

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

1043 [0, b, 0],

1044 [0, 0, 0.]])

1045

1046

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

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

1049class CRECT(BravaisLattice):

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

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

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

1053 # like in ordinary RECT.

1054 #

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

1056 # identical parameters.

1057 if alpha >= 90:

1058 raise UnconventionalLattice(

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

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

1061

1062 def _cell(self, a, alpha):

1063 x = np.cos(alpha * _degrees)

1064 y = np.sin(alpha * _degrees)

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

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

1067 [0, 0, 0.]])

1068

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

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

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

1072 eta = sina2 / sina

1073 cosa = np.cos(alpha * _degrees)

1074 xi = eta * cosa

1075

1076 points = [[0, 0, 0],

1077 [eta, - eta, 0],

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

1079 [0.5, 0.5, 0]]

1080

1081 return points

1082

1083

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

1085 [['SQR', 'GMX', 'MGXM',

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

1087 ndim=2)

1088class SQR(BravaisLattice):

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

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

1091

1092 def _cell(self, a):

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

1094 [0, a, 0],

1095 [0, 0, 0.]])

1096

1097

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

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

1100 ndim=1)

1101class LINE(BravaisLattice):

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

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

1104

1105 def _cell(self, a):

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

1107 [0.0, 0.0, 0.0],

1108 [0.0, 0.0, 0.0]])

1109

1110

1111def celldiff(cell1, cell2):

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

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

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

1115 v1v2 = cell1.volume * cell2.volume

1116 if v1v2 < 1e-10:

1117 # (Proposed cell may be linearly dependent)

1118 return np.inf

1119

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

1121 x1 = cell1 @ cell1.T

1122 x2 = cell2 @ cell2.T

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

1124 return dev

1125

1126

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

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

1129

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

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

1132

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

1134 lattice, raise RuntimeError."""

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

1136

1137

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

1139 """Find Bravais lattice representing this cell.

1140

1141 Returns Bravais lattice object representing the cell along with

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

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

1144 from ase.geometry.bravais_type_engine import niggli_op_table

1145

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

1147 npbc = sum(pbc)

1148

1149 cell = cell.uncomplete(pbc)

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

1151

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

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

1154 # in the table:

1155 memory = {}

1156

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

1158 # the first one we find:

1159 for latname in LatticeChecker.check_orders[npbc]:

1160 # There may be multiple Niggli operations that produce valid

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

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

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

1164 matching_lattices = []

1165

1166 for op_key in niggli_op_table[latname]:

1167 checker_and_op = memory.get(op_key)

1168 if checker_and_op is None:

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

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

1171 checker = LatticeChecker(candidate, eps=eps)

1172 memory[op_key] = (checker, normalization_op)

1173 else:

1174 checker, normalization_op = checker_and_op

1175

1176 lat = checker.query(latname)

1177 if lat is not None:

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

1179 matching_lattices.append((lat, op))

1180

1181 if not matching_lattices:

1182 continue # Move to next Bravais lattice

1183

1184 lat, op = pick_best_lattice(matching_lattices)

1185

1186 if npbc == 2 and op[2, 2] < 0:

1187 op = flip_2d_handedness(op)

1188

1189 return lat, op

1190

1191 raise RuntimeError('Failed to recognize lattice')

1192

1193

1194def flip_2d_handedness(op):

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

1196 # components are necessarily also left-handed which

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

1198 #

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

1200 # z axis and interchanges x/y:

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

1202 return repair_op @ op

1203

1204

1205def pick_best_lattice(matching_lattices):

1206 """Return (lat, op) with lowest orthogonality defect."""

1207 best = None

1208 best_defect = np.inf

1209 for lat, op in matching_lattices:

1210 cell = lat.tocell().complete()

1211 orthogonality_defect = np.prod(cell.lengths()) / cell.volume

1212 if orthogonality_defect < best_defect:

1213 best = lat, op

1214 best_defect = orthogonality_defect

1215 return best

1216

1217

1218class LatticeChecker:

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

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

1221 check_orders = {

1222 1: ['LINE'],

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

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

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

1226

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

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

1229

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

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

1232 by directly through one of the Bravais lattice classes.

1233

1234 Generally for internal use (this module).

1235

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

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

1238 the tolerance eps is not met."""

1239 self.cell = cell

1240 self.eps = eps

1241

1242 self.cellpar = cell.cellpar()

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

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

1245

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

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

1248

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

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

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

1252

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

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

1255 return None

1256 try:

1257 lat = latcls(*args)

1258 except UnconventionalLattice:

1259 return None

1260

1261 newcell = lat.tocell()

1262 err = celldiff(self.cell, newcell)

1263 if err < self.eps:

1264 return lat

1265

1266 def match(self):

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

1268

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

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

1271 lat = self.query(name)

1272 if lat:

1273 return lat

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

1275 'with lengths and angles {}'

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

1277

1278 def query(self, latname):

1279 """Match cell against named Bravais lattice.

1280

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

1282 meth = getattr(self, latname)

1283 lat = meth()

1284 return lat

1285

1286 def LINE(self):

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

1288

1289 def SQR(self):

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

1291

1292 def RECT(self):

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

1294

1295 def CRECT(self):

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

1297

1298 def HEX2D(self):

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

1300

1301 def OBL(self):

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

1303

1304 def CUB(self):

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

1306 # it matches, else None.

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

1308

1309 def FCC(self):

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

1311

1312 def BCC(self):

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

1314

1315 def TET(self):

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

1317

1318 def _bct_orci_lengths(self):

1319 # Coordinate-system independent relation for BCT and ORCI

1320 # standard cells:

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

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

1323 # == b² / 2 (ORCI)

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

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

1326 prods = self.prods

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

1328 if any(lengthsqr < 0):

1329 return None

1330 return np.sqrt(lengthsqr)

1331

1332 def BCT(self):

1333 lengths = self._bct_orci_lengths()

1334 if lengths is None:

1335 return None

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

1337

1338 def HEX(self):

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

1340

1341 def RHL(self):

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

1343

1344 def ORC(self):

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

1346

1347 def ORCF(self):

1348 # ORCF standard cell:

1349 # a2 · a3 = a²/4

1350 # a3 · a1 = b²/4

1351 # a1 · a2 = c²/4

1352 prods = self.prods

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

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

1355 return self._check(ORCF, *orcf_abc)

1356

1357 def ORCI(self):

1358 lengths = self._bct_orci_lengths()

1359 if lengths is None:

1360 return None

1361 return self._check(ORCI, *lengths)

1362

1363 def _orcc_ab(self):

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

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

1366 prods = self.prods

1367 orcc_sqr_ab = np.empty(2)

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

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

1370 if all(orcc_sqr_ab > 0):

1371 return np.sqrt(orcc_sqr_ab)

1372

1373 def ORCC(self):

1374 orcc_lengths_ab = self._orcc_ab()

1375 if orcc_lengths_ab is None:

1376 return None

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

1378

1379 def MCL(self):

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

1381

1382 def MCLC(self):

1383 # MCLC is similar to ORCC:

1384 orcc_ab = self._orcc_ab()

1385 if orcc_ab is None:

1386 return None

1387

1388 prods = self.prods

1389 C = self.C

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

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

1392 if -1 < mclc_cosa < 1:

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

1394 if mclc_b > C:

1395 # XXX Temporary fix for certain otherwise

1396 # unrecognizable lattices.

1397 #

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

1399 # something just outside the domain of conventional

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

1401 # propose a nearby conventional lattice instead, which

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

1403 mclc_b = 0.5 * (mclc_b + C)

1404 C = mclc_b

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

1406

1407 def TRI(self):

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

1409

1410

1411def all_variants():

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

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

1414 alpha = 55.0

1415 yield CUB(a)

1416 yield FCC(a)

1417 yield BCC(a)

1418 yield TET(a, c)

1419 bct1 = BCT(2 * a, c)

1420 bct2 = BCT(a, c)

1421 assert bct1.variant == 'BCT1'

1422 assert bct2.variant == 'BCT2'

1423

1424 yield bct1

1425 yield bct2

1426

1427 yield ORC(a, b, c)

1428

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

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

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

1432 orcf3 = ORCF(a0, b, c)

1433 assert orcf1.variant == 'ORCF1'

1434 assert orcf2.variant == 'ORCF2'

1435 assert orcf3.variant == 'ORCF3'

1436 yield orcf1

1437 yield orcf2

1438 yield orcf3

1439

1440 yield ORCI(a, b, c)

1441 yield ORCC(a, b, c)

1442

1443 yield HEX(a, c)

1444

1445 rhl1 = RHL(a, alpha=55.0)

1446 assert rhl1.variant == 'RHL1'

1447 yield rhl1

1448

1449 rhl2 = RHL(a, alpha=105.0)

1450 assert rhl2.variant == 'RHL2'

1451 yield rhl2

1452

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

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

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

1456

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

1458 assert mclc1.variant == 'MCLC1'

1459 yield mclc1

1460 # mclc2 has same special points as mclc1

1461

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

1463 assert mclc3.variant == 'MCLC3'

1464 yield mclc3

1465 # mclc4 has same special points as mclc3

1466

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

1468

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

1470 assert mclc5.variant == 'MCLC5'

1471 yield mclc5

1472

1473 def get_tri(kcellpar):

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

1475 icell = Cell.fromcellpar(kcellpar)

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

1477 return TRI(*cellpar)

1478

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

1480 assert tri1a.variant == 'TRI1a'

1481 yield tri1a

1482

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

1484 assert tri1b.variant == 'TRI1b'

1485 yield tri1b

1486

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

1488 assert tri2a.variant == 'TRI2a'

1489 yield tri2a

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

1491 assert tri2b.variant == 'TRI2b'

1492 yield tri2b

1493

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

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

1496 # of the same lattice.

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

1498 yield RECT(a, b)

1499 yield CRECT(a, alpha=alpha)

1500 yield HEX2D(a)

1501 yield SQR(a)

1502 yield LINE(a)

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

742 statements