Coverage for icet/tools/geometry.py: 97%
99 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-12-26 04:12 +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) -> List[str]:
241 """Returns the Wyckoff symbols of the input structure. The Wyckoff
242 sites are of general interest for symmetry analysis but can be
243 especially useful when setting up, e.g., a
244 :class:`SiteOccupancyObserver
245 <mchammer.observers.SiteOccupancyObserver>`.
246 The Wyckoff labels can be conveniently attached as an array to the
247 structure object as demonstrated in the examples section below.
249 By default the occupation of the sites is part of the symmetry
250 analysis. If a chemically disordered structure is provided this
251 will usually reduce the symmetry substantially. If one is
252 interested in the symmetry of the underlying structure one can
253 control how occupations are handled. To this end, one can provide
254 the :attr:`map_occupations` keyword argument. The latter must be a
255 list, each entry of which is a list of species that should be
256 treated as indistinguishable. As a shortcut, if *all* species
257 should be treated as indistinguishable one can provide an empty
258 list. Examples that illustrate the usage of the keyword are given
259 below.
261 Parameters
262 ----------
263 structure
264 Input structure. Note that the occupation of the sites is
265 included in the symmetry analysis.
266 map_occupations
267 Each sublist in this list specifies a group of chemical
268 species that shall be treated as indistinguishable for the
269 purpose of the symmetry analysis.
270 symprec
271 Tolerance imposed when analyzing the symmetry using spglib.
273 Examples
274 --------
275 Wyckoff sites of a hexagonal-close packed structure::
277 >>> from ase.build import bulk
278 >>> structure = bulk('Ti')
279 >>> wyckoff_sites = get_wyckoff_sites(structure)
280 >>> print(wyckoff_sites)
281 ['2d', '2d']
284 The Wyckoff labels can also be attached as an array to the
285 structure, in which case the information is also included when
286 storing the Atoms object::
288 >>> from ase.io import write
289 >>> structure.new_array('wyckoff_sites', wyckoff_sites, str)
290 >>> write('structure.xyz', structure)
292 The function can also be applied to supercells::
294 >>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2)
295 >>> wyckoff_sites = get_wyckoff_sites(structure)
296 >>> print(wyckoff_sites)
297 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
298 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
300 Now assume that one is given a supercell of a (Ga,Al)As
301 alloy. Applying the function directly yields much lower symmetry
302 since the symmetry of the original structure is broken::
304 >>> structure.set_chemical_symbols(
305 ... ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As',
306 ... 'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As'])
307 >>> print(get_wyckoff_sites(structure))
308 ['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i',
309 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i']
311 Since Ga and Al occupy the same sublattice, they should, however,
312 be treated as indistinguishable for the purpose of the symmetry
313 analysis, which can be achieved via the :attr:`map_occupations`
314 keyword::
316 >>> print(get_wyckoff_sites(structure, map_occupations=[['Ga', 'Al'], ['As']]))
317 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
318 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
320 If occupations are to ignored entirely, one can simply provide an
321 empty list. In the present case, this turns the zincblende lattice
322 into a diamond lattice, on which case there is only one Wyckoff
323 site::
325 >>> print(get_wyckoff_sites(structure, map_occupations=[]))
326 ['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a',
327 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a']
328 """
329 structure_copy = structure.copy()
330 if map_occupations is not None:
331 if len(map_occupations) > 0:
332 new_symbols = []
333 for symb in structure_copy.get_chemical_symbols():
334 for group in map_occupations: 334 ↛ 333line 334 didn't jump to line 333, because the loop on line 334 didn't complete
335 if symb in group:
336 new_symbols.append(group[0])
337 break
338 else:
339 new_symbols = len(structure) * ['H']
340 structure_copy.set_chemical_symbols(new_symbols)
341 dataset = spglib.get_symmetry_dataset(ase_atoms_to_spglib_cell(structure_copy), symprec=symprec)
342 n_unitcells = np.linalg.det(dataset['transformation_matrix'])
344 equivalent_atoms = list(dataset['equivalent_atoms'])
345 wyckoffs = {}
346 for index in set(equivalent_atoms):
347 multiplicity = list(dataset['equivalent_atoms']).count(index) / n_unitcells
348 multiplicity = int(round(multiplicity))
349 wyckoffs[index] = '{}{}'.format(multiplicity, dataset['wyckoffs'][index])
351 return [wyckoffs[equivalent_atoms[a.index]] for a in structure_copy]