Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1from itertools import product 

2from typing import Dict, List, Tuple, Generator, Union 

3 

4import numpy as np 

5 

6from ase import Atoms 

7from ase.build import make_supercell 

8from spglib import get_symmetry 

9from spglib import niggli_reduce as spg_nigg_red 

10 

11from icet.tools.geometry import ase_atoms_to_spglib_cell 

12from .structure_enumeration_support.normal_form_matrices \ 

13 import HermiteNormalForm, yield_reduced_hnfs, SmithNormalForm, get_unique_snfs 

14from .structure_enumeration_support.labeling_generation \ 

15 import LabelingGenerator 

16 

17 

18def _translate_labelings( 

19 labeling: tuple, 

20 snf: SmithNormalForm, 

21 nsites: int, 

22 include_self: bool = False) -> Generator[Tuple[int], None, None]: 

23 """ 

24 Yields labelings that are equivalent to the original labeling under 

25 translations as dictated by the Smith normal form (SNF) provided. 

26 

27 Parameters 

28 ---------- 

29 labeling 

30 labeling to be translated 

31 snf 

32 SmithNormalForm object 

33 nsites 

34 number of sites in the primitive cell 

35 include_self 

36 if True original labeling will be included 

37 """ 

38 

39 # Compute size of each block within which translations occur 

40 sizes = [nsites * block for block in snf.blocks] 

41 

42 # Loop over all possible translations within group as defined by snf 

43 for trans in product(range(snf.S[0]), range(snf.S[1]), range(snf.S[2])): 

44 if not include_self and sum(trans) == 0: 

45 continue 

46 

47 labeling_trans = () 

48 for i in range(snf.S[0]): 

49 group = (i + trans[0]) % snf.S[0] 

50 block_i = labeling[sizes[0] * group:sizes[0] * (group + 1)] 

51 for j in range(snf.S[1]): 

52 group = (j + trans[1]) % snf.S[1] 

53 block_j = block_i[sizes[1] * group:sizes[1] * (group + 1)] 

54 for k in range(snf.S[2]): 

55 group = (k + trans[2]) % snf.S[2] 

56 labeling_trans += tuple(block_j[sizes[2] * group: 

57 sizes[2] * (group + 1)]) 

58 yield labeling_trans 

59 

60 

61def _get_all_labelings(snf: SmithNormalForm, 

62 labeling_generator: LabelingGenerator, 

63 nsites: int) -> List[tuple]: 

64 """ 

65 Returns inequivalent labelings corresponding to a Smith normal form 

66 (SNF) matrix. Superperiodic labelings as well as labelings that are 

67 equivalent under translations for this particular SNF will not be 

68 included. However, labelings that are equivalent by rotations that 

69 leave the cell (but not the labeling) unchanged will still be 

70 included, since these have to be removed for each Hermite normal 

71 form (HNF) separately. 

72 

73 Parameters 

74 ---------- 

75 snf 

76 SmithNormalForm object 

77 labeling_generator 

78 LabelingGenerator object 

79 nsites 

80 number of sites per primitive cell 

81 """ 

82 labelings = [] 

83 for labeling in labeling_generator.yield_labelings(snf.ncells): 

84 unique = True 

85 for labeling_trans in _translate_labelings(labeling, snf, nsites, 

86 include_self=False): 

87 # Check whether it translates into itself. If so, 

88 # then it has been added with a smaller cell. 

89 if labeling == labeling_trans: 

90 unique = False 

91 break 

92 

93 # Check with previous labelings, 

94 # if labeling can be translated into a previously 

95 # added labeling, then it is not unique 

96 if labeling_trans in labelings: 

97 unique = False 

98 break 

99 if unique: 

100 labelings.append(labeling) 

101 return labelings 

102 

103 

104def _permute_labeling(labeling: tuple, snf: SmithNormalForm, 

105 transformation: List[np.ndarray], 

106 nsites: int) -> Tuple[int]: 

107 """ 

108 Returns permuted labeling according to transformations defined by 

109 transformation. 

110 

111 Parameters 

112 ---------- 

113 labeling 

114 labeling to be rotated 

115 snf 

116 SmithNormalForm object 

117 transformation 

118 transformation consisting of rotation, translation and basis 

119 shift 

120 nsites 

121 number of sites in the primitive cell 

122 """ 

123 

124 # Calculate transformation imposed by LRL multiplication 

125 new_group_order = np.dot(snf.group_order, transformation[0].T) 

126 

127 # Loop over every atom to find its new position 

128 labeling_new = [0] * len(labeling) 

129 for member_index, member in enumerate(new_group_order): 

130 

131 # Transform according to Gp, 

132 # but each site in the primitive cell also transforms in its own way 

133 for basis in range(nsites): 

134 new_cell = member + transformation[1][basis] 

135 

136 # Calculate new index, first by finding the right block, 

137 # then by adding the basis index to that block 

138 new_index = 0 

139 for i in range(3): 

140 new_index += (new_cell[i] % snf.S[i]) * snf.blocks[i] * nsites 

141 new_index += transformation[2][basis] 

142 

143 # Add the contribution to the hash key 

144 element = labeling[member_index * nsites + basis] 

145 labeling_new[new_index] = element 

146 return tuple(labeling_new) 

147 

148 

149def _yield_unique_labelings(labelings: List[int], snf: SmithNormalForm, 

150 hnf: HermiteNormalForm, nsites: int) -> tuple: 

151 """ 

152 Yields labelings that are unique in every imaginable sense. 

153 

154 Parameters 

155 ---------- 

156 labelkeys 

157 list of hash keys to labelings that may still contain labelings 

158 that are equivalent under rotations that leave the supercell 

159 shape unchanged 

160 snf 

161 SmithNormalForm object 

162 hnf 

163 HermiteNormalForm object 

164 nsites 

165 number of sites in the primitive cell 

166 """ 

167 saved_labelings = [] 

168 for labeling in labelings: 

169 

170 # Check whether labeling is just a rotated version of a previous 

171 # labeling. Apply transformation that is specific to the hnf 

172 # and check all translations of the transformed labeling. 

173 unique = True 

174 for transformation in hnf.transformations: 

175 

176 labeling_rot = _permute_labeling(labeling, snf, transformation, 

177 nsites) 

178 

179 # Commonly, the transformation leaves the labeling 

180 # unchanged, so check that first as a special case 

181 # (yields a quite significant speedup) 

182 if labeling_rot == labeling: 

183 continue 

184 

185 # Translate in all possible ways 

186 for labeling_rot_trans in \ 

187 _translate_labelings(labeling_rot, snf, nsites, 

188 include_self=True): 

189 if labeling_rot_trans in saved_labelings: 

190 # Then we have rotated and translated the labeling 

191 # into one that was already yielded 

192 unique = False 

193 break 

194 if not unique: 

195 break 

196 if unique: 

197 # Then we have finally found a unique structure 

198 # defined by an HNF matrix and a labeling 

199 saved_labelings.append(labeling) 

200 yield labeling 

201 

202 

203def _labeling_to_ase_atoms(labeling: tuple, hnf: np.ndarray, cell: np.ndarray, 

204 new_cell: np.ndarray, 

205 basis: np.ndarray, chemical_symbols: List[str], 

206 pbc: List[bool]) -> Atoms: 

207 """ 

208 Returns structure object corresponding to the given labeling using 

209 labeling, HNF matrix and parent lattice. 

210 

211 Parameters 

212 --------- 

213 labeling 

214 permutation of index of elements 

215 hnf 

216 HNF object defining the supercell 

217 cell 

218 basis vectors of primtive cell listed row-wise 

219 new_cell 

220 new cell shape 

221 basis 

222 scaled coordinates to all sites in the primitive cell 

223 chemical_symbols 

224 list of elements, e.g. ``['Au', 'Ag']`` 

225 pbc 

226 periodic boundary conditions of the primitive structure 

227 """ 

228 symbols = [] 

229 positions = [] 

230 count = 0 

231 for i in range(hnf.H[0, 0]): 

232 coord = i * hnf.H[1, 0] 

233 offset10 = coord // hnf.H[0, 0] + coord % hnf.H[0, 0] 

234 coord = i * hnf.H[2, 0] 

235 offset20 = coord // hnf.H[0, 0] + coord % hnf.H[0, 0] 

236 for j in range(hnf.H[1, 1]): 

237 coord = j * hnf.H[2, 1] 

238 offset21 = coord // hnf.H[1, 1] + coord % hnf.H[1, 1] 

239 for k in range(hnf.H[2, 2]): 

240 for basis_vector in basis: 

241 positions.append(i * cell[0] + 

242 (j + offset10) * cell[1] + 

243 (k + offset20 + offset21) * cell[2] + 

244 np.dot(cell.T, basis_vector)) 

245 symbols.append(chemical_symbols[labeling[count]]) 

246 count += 1 

247 structure = Atoms(symbols, positions, cell=new_cell, pbc=(True, True, True)) 

248 structure.wrap() 

249 structure.pbc = pbc 

250 return structure 

251 

252 

253def _get_symmetry_operations(structure: Atoms, 

254 symprec: float, 

255 position_tolerance: float) -> Dict[str, list]: 

256 """ 

257 Returns symmetry operations permissable for a given structure as 

258 obtained via `spglib <https://atztogo.github.io/spglib/>`_. The 

259 symmetry operations consist of three parts: rotation, translation 

260 and basis shifts. The latter define the way that sublattices 

261 shift upon rotation (correponds to `d_Nd` in [HarFor09]_). 

262 

263 Parameters 

264 ---------- 

265 structure 

266 structure for which symmetry operations are sought 

267 symprec 

268 tolerance imposed when analyzing the symmetry using spglib 

269 position_tolerance 

270 tolerance applied when comparing positions in Cartesian coordinates 

271 """ 

272 symmetries = get_symmetry(ase_atoms_to_spglib_cell(structure), symprec=symprec) 

273 assert symmetries, ('spglib.get_symmetry() failed. Please make sure that' 

274 ' the structure object is sensible.') 

275 rotations = symmetries['rotations'] 

276 translations = symmetries['translations'] 

277 

278 basis = structure.get_scaled_positions() 

279 

280 # Calculate how atoms within the primitive cell are shifted (from one site 

281 # to another) and translated (from one primtive cell to another) upon 

282 # operation with rotation matrix. Note that the translations are needed 

283 # here because different sites translate differently. 

284 basis_shifts = np.zeros((len(rotations), len(structure)), dtype='int64') 

285 sites_translations = [] 

286 for i, rotation in enumerate(rotations): 

287 translation = translations[i] 

288 site_translations = [] 

289 for j, basis_element in enumerate(basis): 

290 

291 # Calculate how the site is transformed when operated on by 

292 # symmetry of parent lattice (rotation and translation) 

293 site_rot_trans = np.dot(rotation, basis_element) + translation 

294 

295 # The site may now have been moved to a different site in a 

296 # different cell. We want to separate the two. (In HarFor09, 

297 # basis_shift (site_rot_trans) corresponds to d_Nd and 

298 # site_translation to t_Nd) 

299 site_translation = [0, 0, 0] 

300 for index in range(3): 

301 while site_rot_trans[index] < -position_tolerance: 

302 site_rot_trans[index] += 1 

303 site_translation[index] -= 1 

304 while site_rot_trans[index] > 1 - position_tolerance: 

305 site_rot_trans[index] -= 1 

306 site_translation[index] += 1 

307 site_translations.append(site_translation) 

308 

309 # Find the basis element that the shifted basis correponds to 

310 found = False 

311 for basis_index, basis_element_comp in enumerate(basis): 

312 distance = site_rot_trans - basis_element_comp 

313 

314 # Make sure that they do not differ with a basis vector 

315 for dist_comp_i, dist_comp in enumerate(distance): 

316 if abs(abs(dist_comp) - 1) < position_tolerance: 316 ↛ 317line 316 didn't jump to line 317, because the condition on line 316 was never true

317 distance[dist_comp_i] = 0 

318 

319 if (abs(distance) < position_tolerance).all(): 

320 assert not found 

321 basis_shifts[i, j] = basis_index 

322 found = True 

323 assert found 

324 

325 sites_translations.append(np.array(site_translations)) 

326 

327 symmetries['translations'] = sites_translations 

328 symmetries['basis_shifts'] = basis_shifts 

329 return symmetries 

330 

331 

332def enumerate_structures(structure: Atoms, 

333 sizes: Union[List[int], range], 

334 chemical_symbols: list, 

335 concentration_restrictions: dict = None, 

336 niggli_reduce: bool = None, 

337 symprec: float = 1e-5, 

338 position_tolerance: float = None) -> Atoms: 

339 """ 

340 Yields a sequence of enumerated structures. The function generates 

341 *all* inequivalent structures that are permissible given a certain 

342 lattice. Using the ``chemical_symbols`` and 

343 ``concentration_restrictions`` keyword arguments it is possible to 

344 specify which chemical_symbols are to be included on which site and 

345 in which concentration range. 

346 

347 The function is sensitive to the boundary conditions of the input 

348 structure. An enumeration of, for example, a surface can thus be 

349 performed by setting ``structure.pbc = [True, True, False]``. 

350 

351 The algorithm implemented here was developed by Gus L. W. Hart and 

352 Rodney W. Forcade in Phys. Rev. B **77**, 224115 (2008) 

353 [HarFor08]_ and Phys. Rev. B **80**, 014120 (2009) [HarFor09]_. 

354 

355 Parameters 

356 ---------- 

357 structure 

358 primitive structure from which derivative superstructures should 

359 be generated 

360 sizes 

361 number of sites (included in enumeration) 

362 chemical_symbols 

363 chemical species with which to decorate the structure, e.g., 

364 ``['Au', 'Ag']``; see below for more examples 

365 concentration_restrictions 

366 allowed concentration range for one or more element in 

367 `chemical_symbols`, e.g., ``{'Au': (0, 0.2)}`` will only 

368 enumerate structures in which the Au content is between 0 and 

369 20 %; here, concentration is always defined as the number of 

370 atoms of the specified kind divided by the number of *all* 

371 atoms. 

372 niggli_reduction 

373 if True perform a Niggli reduction with spglib for each 

374 structure; the default is ``True`` if ``structure`` is periodic in 

375 all directions, ``False`` otherwise. 

376 symprec 

377 tolerance imposed when analyzing the symmetry using spglib 

378 position_tolerance 

379 tolerance applied when comparing positions in Cartesian coordinates; 

380 by default this value is set equal to `symprec` 

381 

382 Examples 

383 -------- 

384 

385 The following code snippet illustrates how to enumerate structures 

386 with up to 6 atoms in the unit cell for a binary alloy without any 

387 constraints:: 

388 

389 >>> from ase.build import bulk 

390 >>> prim = bulk('Ag') 

391 >>> for structure in enumerate_structures(structure=prim, 

392 ... sizes=range(1, 5), 

393 ... chemical_symbols=['Ag', 'Au']): 

394 ... pass # Do something with the structure 

395 

396 To limit the concentration range to 10 to 40% Au the code should 

397 be modified as follows:: 

398 

399 >>> conc_restr = {'Au': (0.1, 0.4)} 

400 >>> for structure in enumerate_structures(structure=prim, 

401 ... sizes=range(1, 5), 

402 ... chemical_symbols=['Ag', 'Au'], 

403 ... concentration_restrictions=conc_restr): 

404 ... pass # Do something with the structure 

405 

406 Often one would like to consider mixing on only one 

407 sublattice. This can be achieved as illustrated for a 

408 Ga(1-x)Al(x)As alloy as follows:: 

409 

410 >>> prim = bulk('GaAs', crystalstructure='zincblende', a=5.65) 

411 >>> for structure in enumerate_structures(structure=prim, 

412 ... sizes=range(1, 9), 

413 ... chemical_symbols=[['Ga', 'Al'], ['As']]): 

414 ... pass # Do something with the structure 

415 

416 """ 

417 

418 if position_tolerance is None: 418 ↛ 421line 418 didn't jump to line 421, because the condition on line 418 was never false

419 position_tolerance = symprec 

420 

421 nsites = len(structure) 

422 basis = structure.get_scaled_positions() 

423 

424 # Construct descriptor of where species are allowed to be 

425 if isinstance(chemical_symbols[0], str): 

426 iter_chemical_symbols = [tuple(range(len(chemical_symbols)))] * nsites 

427 elements = chemical_symbols 

428 elif len(chemical_symbols) == nsites: 

429 assert isinstance(chemical_symbols[0][0], str) 

430 elements = [] 

431 for site in chemical_symbols: 

432 for element in site: 

433 if element not in elements: 

434 elements.append(element) 

435 iter_chemical_symbols = [] 

436 for site in chemical_symbols: 

437 iter_chemical_symbols.append(tuple(elements.index(i) 

438 for i in site)) 

439 else: 

440 raise Exception('chemical_symbols needs to be a list of strings ' 

441 'or a list of list of strings.') 

442 

443 # Adapt concentration restrictions to iter_chemical_symbols 

444 if concentration_restrictions: 

445 concentrations = {} 

446 for key, concentration_range in concentration_restrictions.items(): 

447 assert len(concentration_range) == 2, \ 

448 ('Each concentration range' + 

449 ' needs to be specified as (c_low, c_high)') 

450 if key not in elements: 450 ↛ 451line 450 didn't jump to line 451, because the condition on line 450 was never true

451 raise ValueError('{} found in concentration_restrictions but' 

452 ' not in chemical_symbols'.format(key)) 

453 concentrations[elements.index(key)] = concentration_range 

454 else: 

455 concentrations = None 

456 

457 # Construct labeling generator 

458 labeling_generator = LabelingGenerator(iter_chemical_symbols, 

459 concentrations) 

460 

461 # Niggli reduce by default if all directions have 

462 # periodic boundary conditions 

463 if niggli_reduce is None: 463 ↛ 466line 463 didn't jump to line 466, because the condition on line 463 was never false

464 niggli_reduce = (sum(structure.pbc) == 3) 

465 

466 symmetries = _get_symmetry_operations(structure, 

467 symprec=symprec, 

468 position_tolerance=position_tolerance) 

469 

470 # Loop over each cell size 

471 for ncells in sizes: 

472 if ncells == 0: 472 ↛ 473line 472 didn't jump to line 473, because the condition on line 472 was never true

473 continue 

474 

475 hnfs = list(yield_reduced_hnfs(ncells, symmetries, structure.pbc)) 

476 snfs = get_unique_snfs(hnfs) 

477 

478 for snf in snfs: 

479 labelings = _get_all_labelings(snf, labeling_generator, nsites) 

480 for hnf in snf.hnfs: 

481 if niggli_reduce: 

482 new_cell = spg_nigg_red(np.dot(hnf.H.T, structure.cell)) 

483 if new_cell is None: 483 ↛ 484line 483 didn't jump to line 484, because the condition on line 483 was never true

484 new_cell = np.dot(hnf.H.T, structure.cell) 

485 else: 

486 new_cell = np.dot(hnf.H.T, structure.cell) 

487 for labeling in _yield_unique_labelings(labelings, snf, hnf, 

488 nsites): 

489 yield _labeling_to_ase_atoms(labeling, hnf, structure.cell, 

490 new_cell, basis, elements, 

491 structure.pbc) 

492 

493 

494def enumerate_supercells(structure: Atoms, 

495 sizes: Union[List[int], range], 

496 niggli_reduce: bool = None, 

497 symprec: float = 1e-5, 

498 position_tolerance: float = None) -> Atoms: 

499 """ 

500 Yields a sequence of enumerated supercells. The function generates 

501 *all* inequivalent supercells that are permissible given a certain 

502 lattice. Any supercell can be reduced to one of the supercells 

503 generated. 

504 

505 The function is sensitive to the boundary conditions of the input 

506 structure. An enumeration of, for example, a surface can thus be 

507 performed by setting ``structure.pbc = [True, True, False]``. 

508 

509 The algorithm is based on Gus L. W. Hart and 

510 Rodney W. Forcade in Phys. Rev. B **77**, 224115 (2008) 

511 [HarFor08]_ and Phys. Rev. B **80**, 014120 (2009) [HarFor09]_. 

512 

513 Parameters 

514 ---------- 

515 structure 

516 primitive structure from which supercells should be 

517 generated 

518 sizes 

519 number of sites (included in enumeration) 

520 niggli_reduction 

521 if True perform a Niggli reduction with spglib for each 

522 supercell; the default is ``True`` if ``structure`` is periodic in 

523 all directions, ``False`` otherwise. 

524 symprec 

525 tolerance imposed when analyzing the symmetry using spglib 

526 position_tolerance 

527 tolerance applied when comparing positions in Cartesian coordinates; 

528 by default this value is set equal to `symprec` 

529 

530 Examples 

531 -------- 

532 

533 The following code snippet illustrates how to enumerate supercells 

534 with up to 6 atoms in the unit cell:: 

535 

536 >>> from ase.build import bulk 

537 >>> prim = bulk('Ag') 

538 >>> for supercell in enumerate_supercells(structure=prim, sizes=range(1, 7)): 

539 ... pass # Do something with the supercell 

540 

541 """ 

542 

543 if position_tolerance is None: 543 ↛ 548line 543 didn't jump to line 548, because the condition on line 543 was never false

544 position_tolerance = symprec 

545 

546 # Niggli reduce by default if all directions have 

547 # periodic boundary conditions 

548 if niggli_reduce is None: 

549 niggli_reduce = (sum(structure.pbc) == 3) 

550 

551 symmetries = _get_symmetry_operations(structure, 

552 symprec=symprec, 

553 position_tolerance=position_tolerance) 

554 

555 for ncells in sizes: 

556 for hnf in yield_reduced_hnfs(ncells, symmetries, structure.pbc): 

557 supercell = make_supercell(structure, hnf.H.T) 

558 if niggli_reduce: 

559 new_cell = spg_nigg_red(np.dot(hnf.H.T, structure.cell)) 

560 if new_cell is None: # Happens when spglib fails to Niggli reduce 560 ↛ 561line 560 didn't jump to line 561, because the condition on line 560 was never true

561 yield supercell 

562 else: 

563 Pprim = np.dot(new_cell, np.linalg.inv(structure.cell)) 

564 yield make_supercell(structure, Pprim) 

565 else: 

566 yield supercell