Coverage for calorine/tools/structures.py: 100%
79 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-10-08 12:52 +0000
1import numpy as np
2from ase import Atoms
3from ase.filters import FrechetCellFilter
4from ase.optimize import BFGS, LBFGS, FIRE, GPMin
5from ase.optimize.sciopt import SciPyFminBFGS
6from ase.units import GPa
8try:
9 import spglib
10 from spglib._spglib import SpglibCppError
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://docs.ase-lib.org/ase/optimize.html>`_
43 scalar_pressure
44 External pressure in GPa.
45 """
46 if structure.calc is None:
47 raise ValueError('Structure has no attached calculator object')
48 if constant_cell:
49 ucf = structure
50 else:
51 ucf = FrechetCellFilter(
52 structure, constant_volume=constant_volume, scalar_pressure=scalar_pressure * GPa)
53 kwargs['logfile'] = kwargs.get('logfile', None)
54 if minimizer == 'bfgs':
55 dyn = BFGS(ucf, **kwargs)
56 dyn.run(fmax=fmax, steps=steps)
57 elif minimizer == 'lbfgs':
58 dyn = LBFGS(ucf, **kwargs)
59 dyn.run(fmax=fmax, steps=steps)
60 elif minimizer == 'bfgs-scipy':
61 dyn = SciPyFminBFGS(ucf, **kwargs)
62 dyn.run(fmax=fmax, steps=steps)
63 elif minimizer == 'fire':
64 dyn = FIRE(ucf, **kwargs)
65 dyn.run(fmax=fmax, steps=steps)
66 elif minimizer == 'gpmin':
67 dyn = GPMin(ucf, **kwargs)
68 dyn.run(fmax=fmax, steps=steps)
69 else:
70 raise ValueError(f'Unknown minimizer: {minimizer}')
73def get_spacegroup(
74 structure: Atoms,
75 symprec: float = 1e-5,
76 angle_tolerance: float = -1.0,
77 style: str = 'international',
78) -> str:
79 """Returns the space group of a structure using spglib.
80 This is a convenience interface to the :func:`get_spacegroup`
81 function of spglib that works directly with ase Atoms objects.
83 Parameters
84 ----------
85 structure
86 Input atomic structure.
87 symprec
88 Tolerance imposed when analyzing the symmetry.
89 angle_tolerance
90 Tolerance imposed when analyzing angles.
91 style
92 Space group notation to be used. Can be ``'international'`` for the
93 interational tables of crystallography (ITC) style (Hermann-Mauguin
94 and ITC number) or ``'Schoenflies'`` for the Schoenflies notation.
95 """
96 if not spglib_available:
97 raise ImportError('\
98 spglib must be available in order for this function to work, \
99 or you are using a too old version (less than 2.7).') # 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, \
143 or you are using a too old version (less than 2.7).') # pragma: no cover
145 structure_tuple = (
146 structure.get_cell(),
147 structure.get_scaled_positions(),
148 structure.numbers)
149 try:
150 result = spglib.standardize_cell(
151 structure_tuple, to_primitive=to_primitive,
152 no_idealize=no_idealize, symprec=symprec)
153 except SpglibCppError:
154 result = None
155 if result is None:
156 raise ValueError('spglib failed to find the primitive cell, maybe caused by large symprec.')
157 lattice, scaled_positions, numbers = result
158 scaled_positions = [np.round(pos, 12) for pos in scaled_positions]
159 structure_prim = Atoms(scaled_positions=scaled_positions,
160 numbers=numbers, cell=lattice, pbc=True)
161 structure_prim.wrap()
163 return structure_prim
166def get_wyckoff_sites(
167 structure: Atoms,
168 map_occupations: list[list[str]] = None,
169 symprec: float = 1e-5,
170 include_representative_atom_index: bool = False,
171) -> list[str]:
172 """Returns the Wyckoff symbols of the input structure.
173 The Wyckoff labels can be conveniently attached as an array to the
174 structure object as demonstrated in the examples section below.
176 By default the occupation of the sites is part of the symmetry
177 analysis. If a chemically disordered structure is provided this
178 will usually reduce the symmetry substantially. If one is
179 interested in the symmetry of the underlying structure one can
180 control how occupations are handled. To this end, one can provide
181 the :attr:`map_occupations` keyword argument. The latter must be a
182 list, each entry of which is a list of species that should be
183 treated as indistinguishable. As a shortcut, if *all* species
184 should be treated as indistinguishable one can provide an empty
185 list. Examples that illustrate the usage of the keyword are given
186 below.
188 Parameters
189 ----------
190 structure
191 Input structure. Note that the occupation of the sites is
192 included in the symmetry analysis.
193 map_occupations
194 Each sublist in this list specifies a group of chemical
195 species that shall be treated as indistinguishable for the
196 purpose of the symmetry analysis.
197 symprec
198 Tolerance imposed when analyzing the symmetry using spglib.
199 include_representative_atom_index
200 If True the index of the first atom in the structure that is
201 representative of the Wyckoff site is included in the symbol.
202 This is in particular useful in cases when there are multiple
203 Wyckoff sites sites with the same Wyckoff letter.
205 Examples
206 --------
207 Wyckoff sites of a hexagonal-close packed structure::
209 >>> from ase.build import bulk
210 >>> structure = bulk('Ti')
211 >>> wyckoff_sites = get_wyckoff_sites(structure)
212 >>> print(wyckoff_sites)
213 ['2d', '2d']
216 The Wyckoff labels can also be attached as an array to the
217 structure, in which case the information is also included when
218 storing the Atoms object::
220 >>> from ase.io import write
221 >>> structure.new_array('wyckoff_sites', wyckoff_sites, str)
222 >>> write('structure.xyz', structure) # xdoctest: +SKIP
224 The function can also be applied to supercells::
226 >>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2)
227 >>> wyckoff_sites = get_wyckoff_sites(structure)
228 >>> print(wyckoff_sites)
229 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
230 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
232 Now assume that one is given a supercell of a (Ga,Al)As
233 alloy. Applying the function directly yields much lower symmetry
234 since the symmetry of the original structure is broken::
236 >>> structure.set_chemical_symbols(
237 ... ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As',
238 ... 'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As'])
239 >>> print(get_wyckoff_sites(structure))
240 ['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i',
241 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i']
243 Since Ga and Al occupy the same sublattice, they should, however,
244 be treated as indistinguishable for the purpose of the symmetry
245 analysis, which can be achieved via the :attr:`map_occupations`
246 keyword::
248 >>> print(get_wyckoff_sites(structure, map_occupations=[['Ga', 'Al'], ['As']]))
249 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
250 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
252 If occupations are to ignored entirely, one can simply provide an
253 empty list. In the present case, this turns the zincblende lattice
254 into a diamond lattice, on which case there is only one Wyckoff
255 site::
257 >>> print(get_wyckoff_sites(structure, map_occupations=[]))
258 ['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a',
259 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a']
260 """
261 if not spglib_available:
262 raise ImportError('\
263 spglib must be available in order for this function to work, \
264 or you are using a too old version (less than 2.7).') # pragma: no cover
266 structure_copy = structure.copy()
267 if map_occupations is not None:
268 if len(map_occupations) > 0:
269 new_symbols = []
270 for symb in structure_copy.get_chemical_symbols():
271 for group in map_occupations: # pragma: no cover - because of the break
272 if symb in group:
273 new_symbols.append(group[0])
274 break
275 else:
276 new_symbols = len(structure) * ['H']
277 structure_copy.set_chemical_symbols(new_symbols)
278 structure_tuple = (
279 structure_copy.get_cell(),
280 structure_copy.get_scaled_positions(),
281 structure_copy.numbers)
282 dataset = spglib.get_symmetry_dataset(structure_tuple, symprec=symprec)
283 n_unitcells = np.linalg.det(dataset.transformation_matrix)
285 equivalent_atoms = list(dataset.equivalent_atoms)
286 wyckoffs = {}
287 for index in set(equivalent_atoms):
288 multiplicity = list(dataset.equivalent_atoms).count(index) / n_unitcells
289 multiplicity = int(round(multiplicity))
290 wyckoff = '{}{}'.format(multiplicity, dataset.wyckoffs[index])
291 if include_representative_atom_index:
292 wyckoff += f'-{index}'
293 wyckoffs[index] = wyckoff
295 return [wyckoffs[equivalent_atoms[a.index]] for a in structure_copy]