Coverage for calorine/tools/structures.py: 100%
80 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-18 13:01 +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 from spglib._spglib import SpglibCppError
12 spglib_available = True
13except ImportError: # pragma: no cover
14 spglib_available = False
17def relax_structure(structure: Atoms,
18 fmax: float = 0.001,
19 steps: int = 500,
20 minimizer: str = 'bfgs',
21 constant_cell: bool = False,
22 constant_volume: bool = False,
23 scalar_pressure: float = 0.0,
24 **kwargs) -> None:
25 """Relaxes the given structure.
27 Parameters
28 ----------
29 structure
30 Atomic configuration to relax.
31 fmax
32 Stop relaxation if the absolute force for all atoms falls below this value.
33 steps
34 Maximum number of relaxation steps the minimizer is allowed to take.
35 minimizer
36 Minimizer to use; possible values: 'bfgs', 'lbfgs', 'fire', 'gpmin', 'bfgs-scipy'.
37 constant_cell
38 If True do not relax the cell or the volume.
39 constant_volume
40 If True relax the cell shape but keep the volume constant.
41 kwargs
42 Keyword arguments to be handed over to the minimizer; possible arguments can be found
43 in the `ASE documentation <https://wiki.fysik.dtu.dk/ase/ase/optimize.html>`_
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 try:
149 result = spglib.standardize_cell(
150 structure_tuple, to_primitive=to_primitive,
151 no_idealize=no_idealize, symprec=symprec)
152 except SpglibCppError:
153 result = None
154 if result is None:
155 raise ValueError('spglib failed to find the primitive cell, maybe caused by large symprec.')
156 lattice, scaled_positions, numbers = result
157 scaled_positions = [np.round(pos, 12) for pos in scaled_positions]
158 structure_prim = Atoms(scaled_positions=scaled_positions,
159 numbers=numbers, cell=lattice, pbc=True)
160 structure_prim.wrap()
162 return structure_prim
165def get_wyckoff_sites(
166 structure: Atoms,
167 map_occupations: List[List[str]] = None,
168 symprec: float = 1e-5,
169 include_representative_atom_index: bool = False,
170) -> List[str]:
171 """Returns the Wyckoff symbols of the input structure.
172 The Wyckoff labels can be conveniently attached as an array to the
173 structure object as demonstrated in the examples section below.
175 By default the occupation of the sites is part of the symmetry
176 analysis. If a chemically disordered structure is provided this
177 will usually reduce the symmetry substantially. If one is
178 interested in the symmetry of the underlying structure one can
179 control how occupations are handled. To this end, one can provide
180 the :attr:`map_occupations` keyword argument. The latter must be a
181 list, each entry of which is a list of species that should be
182 treated as indistinguishable. As a shortcut, if *all* species
183 should be treated as indistinguishable one can provide an empty
184 list. Examples that illustrate the usage of the keyword are given
185 below.
187 Parameters
188 ----------
189 structure
190 Input structure. Note that the occupation of the sites is
191 included in the symmetry analysis.
192 map_occupations
193 Each sublist in this list specifies a group of chemical
194 species that shall be treated as indistinguishable for the
195 purpose of the symmetry analysis.
196 symprec
197 Tolerance imposed when analyzing the symmetry using spglib.
198 include_representative_atom_index
199 If True the index of the first atom in the structure that is
200 representative of the Wyckoff site is included in the symbol.
201 This is in particular useful in cases when there are multiple
202 Wyckoff sites sites with the same Wyckoff letter.
204 Examples
205 --------
206 Wyckoff sites of a hexagonal-close packed structure::
208 >>> from ase.build import bulk
209 >>> structure = bulk('Ti')
210 >>> wyckoff_sites = get_wyckoff_sites(structure)
211 >>> print(wyckoff_sites)
212 ['2d', '2d']
215 The Wyckoff labels can also be attached as an array to the
216 structure, in which case the information is also included when
217 storing the Atoms object::
219 >>> from ase.io import write
220 >>> structure.new_array('wyckoff_sites', wyckoff_sites, str)
221 >>> write('structure.xyz', structure)
223 The function can also be applied to supercells::
225 >>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2)
226 >>> wyckoff_sites = get_wyckoff_sites(structure)
227 >>> print(wyckoff_sites)
228 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
229 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
231 Now assume that one is given a supercell of a (Ga,Al)As
232 alloy. Applying the function directly yields much lower symmetry
233 since the symmetry of the original structure is broken::
235 >>> structure.set_chemical_symbols(
236 ... ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As',
237 ... 'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As'])
238 >>> print(get_wyckoff_sites(structure))
239 ['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i',
240 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i']
242 Since Ga and Al occupy the same sublattice, they should, however,
243 be treated as indistinguishable for the purpose of the symmetry
244 analysis, which can be achieved via the :attr:`map_occupations`
245 keyword::
247 >>> print(get_wyckoff_sites(structure, map_occupations=[['Ga', 'Al'], ['As']]))
248 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
249 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
251 If occupations are to ignored entirely, one can simply provide an
252 empty list. In the present case, this turns the zincblende lattice
253 into a diamond lattice, on which case there is only one Wyckoff
254 site::
256 >>> print(get_wyckoff_sites(structure, map_occupations=[]))
257 ['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a',
258 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a']
259 """
260 if not spglib_available:
261 raise ImportError('\
262 spglib must be available in order for this function to work.') # pragma: no cover
264 structure_copy = structure.copy()
265 if map_occupations is not None:
266 if len(map_occupations) > 0:
267 new_symbols = []
268 for symb in structure_copy.get_chemical_symbols():
269 for group in map_occupations: # pragma: no cover - because of the break
270 if symb in group:
271 new_symbols.append(group[0])
272 break
273 else:
274 new_symbols = len(structure) * ['H']
275 structure_copy.set_chemical_symbols(new_symbols)
276 structure_tuple = (
277 structure_copy.get_cell(),
278 structure_copy.get_scaled_positions(),
279 structure_copy.numbers)
280 dataset = spglib.get_symmetry_dataset(structure_tuple, symprec=symprec)
281 n_unitcells = np.linalg.det(dataset.transformation_matrix)
283 equivalent_atoms = list(dataset.equivalent_atoms)
284 wyckoffs = {}
285 for index in set(equivalent_atoms):
286 multiplicity = list(dataset.equivalent_atoms).count(index) / n_unitcells
287 multiplicity = int(round(multiplicity))
288 wyckoff = '{}{}'.format(multiplicity, dataset.wyckoffs[index])
289 if include_representative_atom_index:
290 wyckoff += f'-{index}'
291 wyckoffs[index] = wyckoff
293 return [wyckoffs[equivalent_atoms[a.index]] for a in structure_copy]