Coverage for icet/tools/geometry.py: 97%
102 statements
« prev ^ index » next coverage.py v7.5.0, created at 2025-07-13 04:13 +0000
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
9Vector = List[float]
10T = TypeVar('T')
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.
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.
32 """
34 fractional = np.linalg.solve(cell.T, positions.T).T
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
44 return fractional
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.
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 If ``True`` 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)
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()
79 return structure_prim
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.
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])
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)
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)
113 return fractional_positions
116def get_position_from_lattice_site(structure: Atoms, lattice_site: LatticeSite):
117 """
118 Returns the corresponding position from the lattice site.
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)
131def fractional_to_cartesian(structure: Atoms,
132 frac_positions: List[Vector]) -> np.ndarray:
133 """
134 Converts fractional positions into Cartesian positions.
136 Parameters
137 ----------
138 structure
139 Input atomic structure.
140 frac_positions
141 Fractional positions.
142 """
143 return np.dot(frac_positions, structure.cell)
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 not of same size.'
153 f'{len(container)} != {len(permutation)}')
154 if len(set(permutation)) != len(permutation):
155 raise Exception
156 return [container[s] for s in permutation]
159def ase_atoms_to_spglib_cell(structure: Atoms) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
160 """Returns a tuple comprising three components, corresponding to cell
161 metric, atomic positions, and atomic species.
163 Parameters
164 ----------
165 structure
166 Input atomic structure.
167 """
168 return (structure.cell, structure.get_scaled_positions(), structure.get_atomic_numbers())
171def get_occupied_primitive_structure(structure: Atoms,
172 allowed_species: List[List[str]],
173 symprec: float) -> Tuple[Atoms, List[Tuple[str, ...]]]:
174 """Returns an occupied primitive structure with hydrogen on
175 sublattice 1, Helium on sublattice 2 and so on
177 Parameters
178 ----------
179 structure
180 Input structure.
181 allowed_species
182 Chemical symbols that are allowed on each site.
183 symprec
184 Tolerance imposed when analyzing the symmetry using spglib.
185 """
186 if len(structure) != len(allowed_species): 186 ↛ 187line 186 didn't jump to line 187, because the condition on line 186 was never true
187 raise ValueError(
188 'structure and chemical symbols need to be the same size.')
189 sorted_symbols = sorted({tuple(sorted(s)) for s in allowed_species})
191 decorated_primitive = structure.copy()
192 for i, sym in enumerate(allowed_species):
193 sublattice = sorted_symbols.index(tuple(sorted(sym))) + 1
194 decorated_primitive[i].symbol = chemical_symbols[sublattice]
196 decorated_primitive = get_primitive_structure(decorated_primitive, symprec=symprec)
197 decorated_primitive.wrap()
198 primitive_chemical_symbols: List[Tuple[str, ...]] = []
199 for atom in decorated_primitive:
200 sublattice = chemical_symbols.index(atom.symbol)
201 primitive_chemical_symbols.append(sorted_symbols[sublattice - 1])
203 for symbols in allowed_species:
204 if tuple(sorted(symbols)) in primitive_chemical_symbols:
205 index = primitive_chemical_symbols.index(tuple(sorted(symbols)))
206 primitive_chemical_symbols[index] = symbols
207 return decorated_primitive, primitive_chemical_symbols
210def atomic_number_to_chemical_symbol(numbers: List[int]) -> List[str]:
211 """Returns the chemical symbols equivalent to the input atomic
212 numbers.
214 Parameters
215 ----------
216 numbers
217 Atomic numbers.
218 """
220 symbols = [chemical_symbols[number] for number in numbers]
221 return symbols
224def chemical_symbols_to_numbers(symbols: List[str]) -> List[int]:
225 """Returns the atomic numbers equivalent to the input chemical
226 symbols.
228 Parameters
229 ----------
230 symbols
231 Chemical symbols.
232 """
234 numbers = [chemical_symbols.index(symbols) for symbols in symbols]
235 return numbers
238def get_wyckoff_sites(structure: Atoms,
239 map_occupations: List[List[str]] = None,
240 symprec: float = 1e-5,
241 include_representative_atom_index: bool = False) -> List[str]:
242 """Returns the Wyckoff symbols of the input structure. The Wyckoff
243 sites are of general interest for symmetry analysis but can be
244 especially useful when setting up, e.g., a
245 :class:`SiteOccupancyObserver
246 <mchammer.observers.SiteOccupancyObserver>`.
247 The Wyckoff labels can be conveniently attached as an array to the
248 structure object as demonstrated in the examples section below.
250 By default the occupation of the sites is part of the symmetry
251 analysis. If a chemically disordered structure is provided this
252 will usually reduce the symmetry substantially. If one is
253 interested in the symmetry of the underlying structure one can
254 control how occupations are handled. To this end, one can provide
255 the :attr:`map_occupations` keyword argument. The latter must be a
256 list, each entry of which is a list of species that should be
257 treated as indistinguishable. As a shortcut, if *all* species
258 should be treated as indistinguishable one can provide an empty
259 list. Examples that illustrate the usage of the keyword are given
260 below.
262 Parameters
263 ----------
264 structure
265 Input structure. Note that the occupation of the sites is
266 included in the symmetry analysis.
267 map_occupations
268 Each sublist in this list specifies a group of chemical
269 species that shall be treated as indistinguishable for the
270 purpose of the symmetry analysis.
271 symprec
272 Tolerance imposed when analyzing the symmetry using spglib.
273 include_representative_atom_index
274 If True the index of the first atom in the structure that is
275 representative of the Wyckoff site is included in the symbol.
276 This is in particular useful in cases when there are multiple
277 Wyckoff sites sites with the same Wyckoff letter.
279 Examples
280 --------
281 Wyckoff sites of a hexagonal-close packed structure::
283 >>> from ase.build import bulk
284 >>> structure = bulk('Ti')
285 >>> wyckoff_sites = get_wyckoff_sites(structure)
286 >>> print(wyckoff_sites)
287 ['2d', '2d']
290 The Wyckoff labels can also be attached as an array to the
291 structure, in which case the information is also included when
292 storing the Atoms object::
294 >>> from ase.io import write
295 >>> structure.new_array('wyckoff_sites', wyckoff_sites, str)
296 >>> write('structure.xyz', structure)
298 The function can also be applied to supercells::
300 >>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2)
301 >>> wyckoff_sites = get_wyckoff_sites(structure)
302 >>> print(wyckoff_sites)
303 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
304 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
306 Now assume that one is given a supercell of a (Ga,Al)As
307 alloy. Applying the function directly yields much lower symmetry
308 since the symmetry of the original structure is broken::
310 >>> structure.set_chemical_symbols(
311 ... ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As',
312 ... 'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As'])
313 >>> print(get_wyckoff_sites(structure))
314 ['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i',
315 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i']
317 Since Ga and Al occupy the same sublattice, they should, however,
318 be treated as indistinguishable for the purpose of the symmetry
319 analysis, which can be achieved via the :attr:`map_occupations`
320 keyword::
322 >>> print(get_wyckoff_sites(structure, map_occupations=[['Ga', 'Al'], ['As']]))
323 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
324 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
326 If occupations are to ignored entirely, one can simply provide an
327 empty list. In the present case, this turns the zincblende lattice
328 into a diamond lattice, on which case there is only one Wyckoff
329 site::
331 >>> print(get_wyckoff_sites(structure, map_occupations=[]))
332 ['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a',
333 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a']
334 """
335 structure_copy = structure.copy()
336 if map_occupations is not None:
337 if len(map_occupations) > 0:
338 new_symbols = []
339 for symb in structure_copy.get_chemical_symbols():
340 for group in map_occupations: 340 ↛ 339line 340 didn't jump to line 339, because the loop on line 340 didn't complete
341 if symb in group:
342 new_symbols.append(group[0])
343 break
344 else:
345 new_symbols = len(structure) * ['H']
346 structure_copy.set_chemical_symbols(new_symbols)
347 dataset = spglib.get_symmetry_dataset(ase_atoms_to_spglib_cell(structure_copy), symprec=symprec)
348 n_unitcells = np.linalg.det(dataset.transformation_matrix)
350 equivalent_atoms = list(dataset.equivalent_atoms)
351 wyckoffs = {}
352 for index in set(equivalent_atoms):
353 multiplicity = list(dataset.equivalent_atoms).count(index) / n_unitcells
354 multiplicity = int(round(multiplicity))
355 wyckoff = f'{multiplicity}{dataset.wyckoffs[index]}'
356 if include_representative_atom_index:
357 wyckoff += f'-{index}'
358 wyckoffs[index] = wyckoff
360 return [wyckoffs[equivalent_atoms[a.index]] for a in structure_copy]