Coverage for calorine/tools/structures.py: 100%
76 statements
« prev ^ index » next coverage.py v7.8.2, created at 2025-06-24 16:11 +0000
1from typing import List
2import numpy as np
3from ase import Atoms
4from ase.filters import FrechetCellFilter
5from ase.optimize import BFGS, LBFGS, FIRE, GPMin
6from ase.optimize.sciopt import SciPyFminBFGS
7from ase.units import GPa
9try:
10 import spglib
11 spglib_available = True
12except ImportError: # pragma: no cover
13 spglib_available = False
16def relax_structure(structure: Atoms,
17 fmax: float = 0.001,
18 steps: int = 500,
19 minimizer: str = 'bfgs',
20 constant_cell: bool = False,
21 constant_volume: bool = False,
22 scalar_pressure: float = 0.0,
23 **kwargs) -> None:
24 """Relaxes the given structure.
26 Parameters
27 ----------
28 structure
29 Atomic configuration to relax.
30 fmax
31 Stop relaxation if the absolute force for all atoms falls below this value.
32 steps
33 Maximum number of relaxation steps the minimizer is allowed to take.
34 minimizer
35 Minimizer to use; possible values: 'bfgs', 'lbfgs', 'fire', 'gpmin', 'bfgs-scipy'.
36 constant_cell
37 If True do not relax the cell or the volume.
38 constant_volume
39 If True relax the cell shape but keep the volume constant.
40 kwargs
41 Keyword arguments to be handed over to the minimizer; possible arguments can be found
42 in the `ASE documentation <https://wiki.fysik.dtu.dk/ase/ase/optimize.html>`_
43 https://wiki.fysik.dtu.dk/ase/ase/filters.html#the-frechetcellfilter-class.
44 scalar_pressure
45 External pressure in GPa.
46 """
47 if structure.calc is None:
48 raise ValueError('Structure has no attached calculator object')
49 if constant_cell:
50 ucf = structure
51 else:
52 ucf = FrechetCellFilter(
53 structure, constant_volume=constant_volume, scalar_pressure=scalar_pressure * GPa)
54 kwargs['logfile'] = kwargs.get('logfile', None)
55 if minimizer == 'bfgs':
56 dyn = BFGS(ucf, **kwargs)
57 dyn.run(fmax=fmax, steps=steps)
58 elif minimizer == 'lbfgs':
59 dyn = LBFGS(ucf, **kwargs)
60 dyn.run(fmax=fmax, steps=steps)
61 elif minimizer == 'bfgs-scipy':
62 dyn = SciPyFminBFGS(ucf, **kwargs)
63 dyn.run(fmax=fmax, steps=steps)
64 elif minimizer == 'fire':
65 dyn = FIRE(ucf, **kwargs)
66 dyn.run(fmax=fmax, steps=steps)
67 elif minimizer == 'gpmin':
68 dyn = GPMin(ucf, **kwargs)
69 dyn.run(fmax=fmax, steps=steps)
70 else:
71 raise ValueError(f'Unknown minimizer: {minimizer}')
74def get_spacegroup(
75 structure: Atoms,
76 symprec: float = 1e-5,
77 angle_tolerance: float = -1.0,
78 style: str = 'international',
79) -> str:
80 """Returns the space group of a structure using spglib.
81 This is a convenience interface to the :func:`get_spacegroup`
82 function of spglib that works directly with ase Atoms objects.
84 Parameters
85 ----------
86 structure
87 Input atomic structure.
88 symprec
89 Tolerance imposed when analyzing the symmetry.
90 angle_tolerance
91 Tolerance imposed when analyzing angles.
92 style
93 Space group notation to be used. Can be ``'international'`` for the
94 interational tables of crystallography (ITC) style (Hermann-Mauguin
95 and ITC number) or ``'Schoenflies'`` for the Schoenflies notation.
96 """
97 if not spglib_available:
98 raise ImportError('\
99 spglib must be available in order for this function to work.') # pragma: no cover
101 if style == 'international':
102 symbol_type = 0
103 elif style == 'Schoenflies':
104 symbol_type = 1
105 else:
106 raise ValueError(f'Unknown value for style: {style}')
108 structure_tuple = (
109 structure.get_cell(),
110 structure.get_scaled_positions(),
111 structure.numbers)
112 spg = spglib.get_spacegroup(
113 structure_tuple, symprec=symprec,
114 angle_tolerance=angle_tolerance, symbol_type=symbol_type)
116 return spg
119def get_primitive_structure(
120 structure: Atoms,
121 no_idealize: bool = True,
122 to_primitive: bool = True,
123 symprec: float = 1e-5,
124) -> Atoms:
125 """Returns the primitive structure using spglib.
126 This is a convenience interface to the :func:`standardize_cell`
127 function of spglib that works directly with ase Atoms objects.
129 Parameters
130 ----------
131 structure
132 Input atomic structure.
133 no_idealize
134 If ``True`` lengths and angles are not idealized.
135 to_primitive
136 If ``True`` convert to primitive structure.
137 symprec
138 Tolerance imposed when analyzing the symmetry.
139 """
140 if not spglib_available:
141 raise ImportError('\
142 spglib must be available in order for this function to work.') # pragma: no cover
144 structure_tuple = (
145 structure.get_cell(),
146 structure.get_scaled_positions(),
147 structure.numbers)
148 result = spglib.standardize_cell(
149 structure_tuple, to_primitive=to_primitive,
150 no_idealize=no_idealize, symprec=symprec)
151 if result is None:
152 raise ValueError('spglib failed to find the primitive cell, maybe caused by large symprec.')
153 lattice, scaled_positions, numbers = result
154 scaled_positions = [np.round(pos, 12) for pos in scaled_positions]
155 structure_prim = Atoms(scaled_positions=scaled_positions,
156 numbers=numbers, cell=lattice, pbc=True)
157 structure_prim.wrap()
159 return structure_prim
162def get_wyckoff_sites(
163 structure: Atoms,
164 map_occupations: List[List[str]] = None,
165 symprec: float = 1e-5,
166 include_representative_atom_index: bool = False,
167) -> List[str]:
168 """Returns the Wyckoff symbols of the input structure.
169 The Wyckoff labels can be conveniently attached as an array to the
170 structure object as demonstrated in the examples section below.
172 By default the occupation of the sites is part of the symmetry
173 analysis. If a chemically disordered structure is provided this
174 will usually reduce the symmetry substantially. If one is
175 interested in the symmetry of the underlying structure one can
176 control how occupations are handled. To this end, one can provide
177 the :attr:`map_occupations` keyword argument. The latter must be a
178 list, each entry of which is a list of species that should be
179 treated as indistinguishable. As a shortcut, if *all* species
180 should be treated as indistinguishable one can provide an empty
181 list. Examples that illustrate the usage of the keyword are given
182 below.
184 Parameters
185 ----------
186 structure
187 Input structure. Note that the occupation of the sites is
188 included in the symmetry analysis.
189 map_occupations
190 Each sublist in this list specifies a group of chemical
191 species that shall be treated as indistinguishable for the
192 purpose of the symmetry analysis.
193 symprec
194 Tolerance imposed when analyzing the symmetry using spglib.
195 include_representative_atom_index
196 If True the index of the first atom in the structure that is
197 representative of the Wyckoff site is included in the symbol.
198 This is in particular useful in cases when there are multiple
199 Wyckoff sites sites with the same Wyckoff letter.
201 Examples
202 --------
203 Wyckoff sites of a hexagonal-close packed structure::
205 >>> from ase.build import bulk
206 >>> structure = bulk('Ti')
207 >>> wyckoff_sites = get_wyckoff_sites(structure)
208 >>> print(wyckoff_sites)
209 ['2d', '2d']
212 The Wyckoff labels can also be attached as an array to the
213 structure, in which case the information is also included when
214 storing the Atoms object::
216 >>> from ase.io import write
217 >>> structure.new_array('wyckoff_sites', wyckoff_sites, str)
218 >>> write('structure.xyz', structure)
220 The function can also be applied to supercells::
222 >>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2)
223 >>> wyckoff_sites = get_wyckoff_sites(structure)
224 >>> print(wyckoff_sites)
225 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
226 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
228 Now assume that one is given a supercell of a (Ga,Al)As
229 alloy. Applying the function directly yields much lower symmetry
230 since the symmetry of the original structure is broken::
232 >>> structure.set_chemical_symbols(
233 ... ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As',
234 ... 'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As'])
235 >>> print(get_wyckoff_sites(structure))
236 ['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i',
237 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i']
239 Since Ga and Al occupy the same sublattice, they should, however,
240 be treated as indistinguishable for the purpose of the symmetry
241 analysis, which can be achieved via the :attr:`map_occupations`
242 keyword::
244 >>> print(get_wyckoff_sites(structure, map_occupations=[['Ga', 'Al'], ['As']]))
245 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
246 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
248 If occupations are to ignored entirely, one can simply provide an
249 empty list. In the present case, this turns the zincblende lattice
250 into a diamond lattice, on which case there is only one Wyckoff
251 site::
253 >>> print(get_wyckoff_sites(structure, map_occupations=[]))
254 ['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a',
255 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a']
256 """
257 if not spglib_available:
258 raise ImportError('\
259 spglib must be available in order for this function to work.') # pragma: no cover
261 structure_copy = structure.copy()
262 if map_occupations is not None:
263 if len(map_occupations) > 0:
264 new_symbols = []
265 for symb in structure_copy.get_chemical_symbols():
266 for group in map_occupations: # pragma: no cover - because of the break
267 if symb in group:
268 new_symbols.append(group[0])
269 break
270 else:
271 new_symbols = len(structure) * ['H']
272 structure_copy.set_chemical_symbols(new_symbols)
273 structure_tuple = (
274 structure_copy.get_cell(),
275 structure_copy.get_scaled_positions(),
276 structure_copy.numbers)
277 dataset = spglib.get_symmetry_dataset(structure_tuple, symprec=symprec)
278 n_unitcells = np.linalg.det(dataset.transformation_matrix)
280 equivalent_atoms = list(dataset.equivalent_atoms)
281 wyckoffs = {}
282 for index in set(equivalent_atoms):
283 multiplicity = list(dataset.equivalent_atoms).count(index) / n_unitcells
284 multiplicity = int(round(multiplicity))
285 wyckoff = '{}{}'.format(multiplicity, dataset.wyckoffs[index])
286 if include_representative_atom_index:
287 wyckoff += f'-{index}'
288 wyckoffs[index] = wyckoff
290 return [wyckoffs[equivalent_atoms[a.index]] for a in structure_copy]