Coverage for /builds/debichem-team/python-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 if not matching_lattices:

1180 continue # Move to next Bravais lattice

1181

1182 lat, op = pick_best_lattice(matching_lattices)

1183

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

1185 op = flip_2d_handedness(op)

1186

1187 return lat, op

1188

1189 raise RuntimeError('Failed to recognize lattice')

1190

1191

1192def flip_2d_handedness(op):

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

1194 # components are necessarily also left-handed which

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

1196 #

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

1198 # z axis and interchanges x/y:

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

1200 return repair_op @ op

1201

1202

1203def pick_best_lattice(matching_lattices):

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

1205 best = None

1206 best_defect = np.inf

1207 for lat, op in matching_lattices:

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

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

1210 if orthogonality_defect < best_defect:

1211 best = lat, op

1212 best_defect = orthogonality_defect

1213 return best

1214

1215

1216class LatticeChecker:

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

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

1219 check_orders = {

1220 1: ['LINE'],

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

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

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

1224

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

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

1227

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

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

1230 by directly through one of the Bravais lattice classes.

1231

1232 Generally for internal use (this module).

1233

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

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

1236 the tolerance eps is not met."""

1237 self.cell = cell

1238 self.eps = eps

1239

1240 self.cellpar = cell.cellpar()

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

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

1243

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

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

1246

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

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

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

1250

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

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

1253 return None

1254 try:

1255 lat = latcls(*args)

1256 except UnconventionalLattice:

1257 return None

1258

1259 newcell = lat.tocell()

1260 err = celldiff(self.cell, newcell)

1261 if err < self.eps:

1262 return lat

1263

1264 def match(self):

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

1266

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

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

1269 lat = self.query(name)

1270 if lat:

1271 return lat

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

1273 'with lengths and angles {}'

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

1275

1276 def query(self, latname):

1277 """Match cell against named Bravais lattice.

1278

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

1280 meth = getattr(self, latname)

1281 lat = meth()

1282 return lat

1283

1284 def LINE(self):

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

1286

1287 def SQR(self):

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

1289

1290 def RECT(self):

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

1292

1293 def CRECT(self):

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

1295

1296 def HEX2D(self):

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

1298

1299 def OBL(self):

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

1301

1302 def CUB(self):

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

1304 # it matches, else None.

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

1306

1307 def FCC(self):

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

1309

1310 def BCC(self):

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

1312

1313 def TET(self):

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

1315

1316 def _bct_orci_lengths(self):

1317 # Coordinate-system independent relation for BCT and ORCI

1318 # standard cells:

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

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

1321 # == b² / 2 (ORCI)

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

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

1324 prods = self.prods

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

1326 if any(lengthsqr < 0):

1327 return None

1328 return np.sqrt(lengthsqr)

1329

1330 def BCT(self):

1331 lengths = self._bct_orci_lengths()

1332 if lengths is None:

1333 return None

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

1335

1336 def HEX(self):

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

1338

1339 def RHL(self):

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

1341

1342 def ORC(self):

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

1344

1345 def ORCF(self):

1346 # ORCF standard cell:

1347 # a2 · a3 = a²/4

1348 # a3 · a1 = b²/4

1349 # a1 · a2 = c²/4

1350 prods = self.prods

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

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

1353 return self._check(ORCF, *orcf_abc)

1354

1355 def ORCI(self):

1356 lengths = self._bct_orci_lengths()

1357 if lengths is None:

1358 return None

1359 return self._check(ORCI, *lengths)

1360

1361 def _orcc_ab(self):

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

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

1364 prods = self.prods

1365 orcc_sqr_ab = np.empty(2)

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

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

1368 if all(orcc_sqr_ab > 0):

1369 return np.sqrt(orcc_sqr_ab)

1370

1371 def ORCC(self):

1372 orcc_lengths_ab = self._orcc_ab()

1373 if orcc_lengths_ab is None:

1374 return None

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

1376

1377 def MCL(self):

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

1379

1380 def MCLC(self):

1381 # MCLC is similar to ORCC:

1382 orcc_ab = self._orcc_ab()

1383 if orcc_ab is None:

1384 return None

1385

1386 prods = self.prods

1387 C = self.C

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

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

1390 if -1 < mclc_cosa < 1:

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

1392 if mclc_b > C:

1393 # XXX Temporary fix for certain otherwise

1394 # unrecognizable lattices.

1395 #

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

1397 # something just outside the domain of conventional

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

1399 # propose a nearby conventional lattice instead, which

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

1401 mclc_b = 0.5 * (mclc_b + C)

1402 C = mclc_b

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

1404

1405 def TRI(self):

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

1407

1408

1409def all_variants():

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

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

1412 alpha = 55.0

1413 yield CUB(a)

1414 yield FCC(a)

1415 yield BCC(a)

1416 yield TET(a, c)

1417 bct1 = BCT(2 * a, c)

1418 bct2 = BCT(a, c)

1419 assert bct1.variant == 'BCT1'

1420 assert bct2.variant == 'BCT2'

1421

1422 yield bct1

1423 yield bct2

1424

1425 yield ORC(a, b, c)

1426

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

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

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

1430 orcf3 = ORCF(a0, b, c)

1431 assert orcf1.variant == 'ORCF1'

1432 assert orcf2.variant == 'ORCF2'

1433 assert orcf3.variant == 'ORCF3'

1434 yield orcf1

1435 yield orcf2

1436 yield orcf3

1437

1438 yield ORCI(a, b, c)

1439 yield ORCC(a, b, c)

1440

1441 yield HEX(a, c)

1442

1443 rhl1 = RHL(a, alpha=55.0)

1444 assert rhl1.variant == 'RHL1'

1445 yield rhl1

1446

1447 rhl2 = RHL(a, alpha=105.0)

1448 assert rhl2.variant == 'RHL2'

1449 yield rhl2

1450

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

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

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

1454

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

1456 assert mclc1.variant == 'MCLC1'

1457 yield mclc1

1458 # mclc2 has same special points as mclc1

1459

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

1461 assert mclc3.variant == 'MCLC3'

1462 yield mclc3

1463 # mclc4 has same special points as mclc3

1464

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

1466

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

1468 assert mclc5.variant == 'MCLC5'

1469 yield mclc5

1470

1471 def get_tri(kcellpar):

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

1473 icell = Cell.fromcellpar(kcellpar)

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

1475 return TRI(*cellpar)

1476

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

1478 assert tri1a.variant == 'TRI1a'

1479 yield tri1a

1480

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

1482 assert tri1b.variant == 'TRI1b'

1483 yield tri1b

1484

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

1486 assert tri2a.variant == 'TRI2a'

1487 yield tri2a

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

1489 assert tri2b.variant == 'TRI2b'

1490 yield tri2b

1491

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

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

1494 # of the same lattice.

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

1496 yield RECT(a, b)

1497 yield CRECT(a, alpha=alpha)

1498 yield HEX2D(a)

1499 yield SQR(a)

1500 yield LINE(a)

Coverage for /builds/debichem-team/python-ase/ase/lattice/init.py: 97.30%

742 statements