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

1import numpy as np 

2import spglib 

3from typing import List, Sequence, Tuple, TypeVar 

4from ase import Atoms 

5from ase.data import chemical_symbols 

6from icet.core.lattice_site import LatticeSite 

7from icet.core.structure import Structure 

8 

9Vector = List[float] 

10T = TypeVar('T') 

11 

12 

13def get_scaled_positions(positions: np.ndarray, 

14 cell: np.ndarray, 

15 wrap: bool = True, 

16 pbc: List[bool] = [True, True, True]) -> np.ndarray: 

17 """ 

18 Returns the positions in reduced (or scaled) coordinates. 

19 

20 Parameters 

21 ---------- 

22 positions 

23 atomic positions in cartesian coordinates 

24 cell 

25 cell metric 

26 wrap 

27 if True, positions outside the unit cell will be wrapped into 

28 the cell in the directions with periodic boundary conditions 

29 such that the scaled coordinates are between zero and one. 

30 pbc 

31 periodic boundary conditions flags 

32 """ 

33 

34 fractional = np.linalg.solve(cell.T, positions.T).T 

35 

36 if wrap: 

37 for i, periodic in enumerate(pbc): 

38 if periodic: 

39 # Yes, we need to do it twice. 

40 # See the scaled_positions.py test. 

41 fractional[:, i] %= 1.0 

42 fractional[:, i] %= 1.0 

43 

44 return fractional 

45 

46 

47def get_primitive_structure(structure: Atoms, 

48 no_idealize: bool = True, 

49 to_primitive: bool = True, 

50 symprec: float = 1e-5) -> Atoms: 

51 """ 

52 Returns the primitive structure using spglib. 

53 

54 Parameters 

55 ---------- 

56 structure 

57 input atomic structure 

58 no_idealize 

59 if True lengths and angles are not idealized 

60 to_primitive 

61 convert to primitive structure 

62 symprec 

63 tolerance imposed when analyzing the symmetry using spglib 

64 """ 

65 structure_copy = structure.copy() 

66 structure_as_tuple = ase_atoms_to_spglib_cell(structure_copy) 

67 

68 ret = spglib.standardize_cell( 

69 structure_as_tuple, to_primitive=to_primitive, 

70 no_idealize=no_idealize, symprec=symprec) 

71 if ret is None: 

72 raise ValueError('spglib failed to find the primitive cell, maybe caused by large symprec.') 

73 lattice, scaled_positions, numbers = ret 

74 scaled_positions = [np.round(pos, 12) for pos in scaled_positions] 

75 structure_prim = Atoms(scaled_positions=scaled_positions, 

76 numbers=numbers, cell=lattice, pbc=structure.pbc) 

77 structure_prim.wrap() 

78 

79 return structure_prim 

80 

81 

82def get_fractional_positions_from_neighbor_list( 

83 structure: Structure, neighbor_list: List) -> List[Vector]: 

84 """ 

85 Returns the fractional positions of the lattice sites in structure from 

86 a neighbor list. 

87 

88 Parameters 

89 ---------- 

90 structure 

91 input atomic structure 

92 neighbor_list 

93 list of lattice neighbors of the input structure 

94 """ 

95 neighbor_positions = [] 

96 fractional_positions = [] 

97 lattice_site = LatticeSite(0, [0, 0, 0]) 

98 

99 for i in range(len(neighbor_list)): 

100 lattice_site.index = i 

101 position = structure.get_position(lattice_site) 

102 neighbor_positions.append(position) 

103 for neighbor in neighbor_list[i]: 

104 position = structure.get_position(neighbor) 

105 neighbor_positions.append(position) 

106 

107 if len(neighbor_positions) > 0: 107 ↛ 113line 107 didn't jump to line 113, because the condition on line 107 was never false

108 fractional_positions = get_scaled_positions( 

109 np.array(neighbor_positions), 

110 structure.cell, wrap=False, 

111 pbc=structure.pbc) 

112 

113 return fractional_positions 

114 

115 

116def get_position_from_lattice_site(structure: Atoms, lattice_site: LatticeSite): 

117 """ 

118 Gets the corresponding position from the lattice site. 

119 

120 Parameters 

121 --------- 

122 structure 

123 input atomic structure 

124 lattice_site 

125 specific lattice site of the input structure 

126 """ 

127 return structure[lattice_site.index].position + \ 

128 np.dot(lattice_site.unitcell_offset, structure.cell) 

129 

130 

131def fractional_to_cartesian(structure: Atoms, 

132 frac_positions: List[Vector]) -> np.ndarray: 

133 """ 

134 Turns fractional positions into cartesian positions. 

135 

136 Parameters 

137 ---------- 

138 structure 

139 input atomic structure 

140 frac_positions 

141 fractional positions 

142 """ 

143 return np.dot(frac_positions, structure.cell) 

144 

145 

146def get_permutation(container: Sequence[T], 

147 permutation: List[int]) -> Sequence[T]: 

148 """ 

149 Returns the permuted version of container. 

150 """ 

151 if len(permutation) != len(container): 

152 raise RuntimeError('Container and permutation' 

153 ' not of same size {} != {}'.format( 

154 len(container), len(permutation))) 

155 if len(set(permutation)) != len(permutation): 

156 raise Exception 

157 return [container[s] for s in permutation] 

158 

159 

160def ase_atoms_to_spglib_cell(structure: Atoms) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: 

161 """ 

162 Returns a tuple of three components: cell metric, atomic positions, and 

163 atomic species of the input ASE Atoms object. 

164 """ 

165 return (structure.cell, structure.get_scaled_positions(), structure.get_atomic_numbers()) 

166 

167 

168def get_occupied_primitive_structure(structure: Atoms, 

169 allowed_species: List[List[str]], 

170 symprec: float) -> Tuple[Atoms, List[Tuple[str, ...]]]: 

171 """ 

172 Returns an occupied primitive structure. 

173 Will put hydrogen on sublattice 1, Helium on sublattice 2 and 

174 so on 

175 

176 Parameters 

177 ---------- 

178 structure 

179 input structure 

180 allowed_species 

181 chemical symbols that are allowed on each site 

182 symprec 

183 tolerance imposed when analyzing the symmetry using spglib 

184 

185 Todo 

186 ---- 

187 simplify the revert back to unsorted symbols 

188 """ 

189 if len(structure) != len(allowed_species): 189 ↛ 190line 189 didn't jump to line 190, because the condition on line 189 was never true

190 raise ValueError( 

191 'structure and chemical symbols need to be the same size.') 

192 sorted_symbols = sorted({tuple(sorted(s)) for s in allowed_species}) 

193 

194 decorated_primitive = structure.copy() 

195 for i, sym in enumerate(allowed_species): 

196 sublattice = sorted_symbols.index(tuple(sorted(sym))) + 1 

197 decorated_primitive[i].symbol = chemical_symbols[sublattice] 

198 

199 decorated_primitive = get_primitive_structure(decorated_primitive, symprec=symprec) 

200 decorated_primitive.wrap() 

201 primitive_chemical_symbols: List[Tuple[str, ...]] = [] 

202 for atom in decorated_primitive: 

203 sublattice = chemical_symbols.index(atom.symbol) 

204 primitive_chemical_symbols.append(sorted_symbols[sublattice - 1]) 

205 

206 for symbols in allowed_species: 

207 if tuple(sorted(symbols)) in primitive_chemical_symbols: 

208 index = primitive_chemical_symbols.index(tuple(sorted(symbols))) 

209 primitive_chemical_symbols[index] = symbols 

210 return decorated_primitive, primitive_chemical_symbols 

211 

212 

213def atomic_number_to_chemical_symbol(numbers: List[int]) -> List[str]: 

214 """Returns the chemical symbols equivalent to the input atomic 

215 numbers. 

216 

217 Parameters 

218 ---------- 

219 numbers 

220 atomic numbers 

221 """ 

222 

223 symbols = [chemical_symbols[number] for number in numbers] 

224 return symbols 

225 

226 

227def chemical_symbols_to_numbers(symbols: List[str]) -> List[int]: 

228 """Returns the atomic numbers equivalent to the input chemical 

229 symbols. 

230 

231 Parameters 

232 ---------- 

233 symbols 

234 chemical symbols 

235 """ 

236 

237 numbers = [chemical_symbols.index(symbols) for symbols in symbols] 

238 return numbers 

239 

240 

241def get_wyckoff_sites(structure: Atoms, 

242 map_occupations: List[List[str]] = None, 

243 symprec: float = 1e-5) -> List[str]: 

244 """Returns the Wyckoff symbols of the input structure. The Wyckoff 

245 sites are of general interest for symmetry analysis but can be 

246 especially useful when setting up, e.g., a 

247 :class:`SiteOccupancyObserver 

248 <mchammer.observers.SiteOccupancyObserver>`. 

249 The Wyckoff labels can be conveniently attached as an array to the 

250 structure object as demonstrated in the examples section below. 

251 

252 By default the occupation of the sites is part of the symmetry 

253 analysis. If a chemically disordered structure is provided this 

254 will usually reduce the symmetry substantially. If one is 

255 interested in the symmetry of the underlying structure one can 

256 control how occupations are handled. To this end, one can provide 

257 the ``map_occupations`` keyword argument. The latter must be a 

258 list, each entry of which is a list of species that should be 

259 treated as indistinguishable. As a shortcut, if *all* species 

260 should be treated as indistinguishable one can provide an empty 

261 list. Examples that illustrate the usage of the keyword are given 

262 below. 

263 

264 Parameters 

265 ---------- 

266 structure 

267 input structure, note that the occupation of the sites is 

268 included in the symmetry analysis 

269 map_occupations 

270 each sublist in this list specifies a group of chemical 

271 species that shall be treated as indistinguishable for the 

272 purpose of the symmetry analysis 

273 symprec 

274 tolerance imposed when analyzing the symmetry using spglib 

275 

276 Examples 

277 -------- 

278 Wyckoff sites of a hexagonal-close packed structure:: 

279 

280 >>> from ase.build import bulk 

281 >>> structure = bulk('Ti') 

282 >>> wyckoff_sites = get_wyckoff_sites(structure) 

283 >>> print(wyckoff_sites) 

284 ['2d', '2d'] 

285 

286 

287 The Wyckoff labels can also be attached as an array to the 

288 structure, in which case the information is also included when 

289 storing the Atoms object:: 

290 

291 >>> from ase.io import write 

292 >>> structure.new_array('wyckoff_sites', wyckoff_sites, str) 

293 >>> write('structure.xyz', structure) 

294 

295 The function can also be applied to supercells:: 

296 

297 >>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2) 

298 >>> wyckoff_sites = get_wyckoff_sites(structure) 

299 >>> print(wyckoff_sites) 

300 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c', 

301 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c'] 

302 

303 Now assume that one is given a supercell of a (Ga,Al)As 

304 alloy. Applying the function directly yields much lower symmetry 

305 since the symmetry of the original structure is broken:: 

306 

307 >>> structure.set_chemical_symbols( 

308 ... ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As', 

309 ... 'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As']) 

310 >>> print(get_wyckoff_sites(structure)) 

311 ['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i', 

312 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i'] 

313 

314 Since Ga and Al occupy the same sublattice, they should, however, 

315 be treated as indistinguishable for the purpose of the symmetry 

316 analysis, which can be achieved via the ``map_occupations`` 

317 keyword:: 

318 

319 >>> print(get_wyckoff_sites(structure, map_occupations=[['Ga', 'Al'], ['As']])) 

320 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c', 

321 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c'] 

322 

323 If occupations are to ignored entirely, one can simply provide an 

324 empty list. In the present case, this turns the zincblende lattice 

325 into a diamond lattice, on which case there is only one Wyckoff 

326 site:: 

327 

328 >>> print(get_wyckoff_sites(structure, map_occupations=[])) 

329 ['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a', 

330 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a'] 

331 """ 

332 structure_copy = structure.copy() 

333 if map_occupations is not None: 

334 if len(map_occupations) > 0: 

335 new_symbols = [] 

336 for symb in structure_copy.get_chemical_symbols(): 

337 for group in map_occupations: 337 ↛ 336line 337 didn't jump to line 336, because the loop on line 337 didn't complete

338 if symb in group: 

339 new_symbols.append(group[0]) 

340 break 

341 else: 

342 new_symbols = len(structure) * ['H'] 

343 structure_copy.set_chemical_symbols(new_symbols) 

344 dataset = spglib.get_symmetry_dataset(ase_atoms_to_spglib_cell(structure_copy), symprec=symprec) 

345 n_unitcells = np.linalg.det(dataset['transformation_matrix']) 

346 

347 equivalent_atoms = list(dataset['equivalent_atoms']) 

348 wyckoffs = {} 

349 for index in set(equivalent_atoms): 

350 multiplicity = list(dataset['equivalent_atoms']).count(index) / n_unitcells 

351 multiplicity = int(round(multiplicity)) 

352 wyckoffs[index] = '{}{}'.format(multiplicity, dataset['wyckoffs'][index]) 

353 

354 return [wyckoffs[equivalent_atoms[a.index]] for a in structure_copy]