Coverage for calorine/tools/structures.py: 100%
80 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-23 12:50 +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://docs.ase-lib.org/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, \
100 or you are using a too old version (less than 2.7).') # pragma: no cover
102 if style == 'international':
103 symbol_type = 0
104 elif style == 'Schoenflies':
105 symbol_type = 1
106 else:
107 raise ValueError(f'Unknown value for style: {style}')
109 structure_tuple = (
110 structure.get_cell(),
111 structure.get_scaled_positions(),
112 structure.numbers)
113 spg = spglib.get_spacegroup(
114 structure_tuple, symprec=symprec,
115 angle_tolerance=angle_tolerance, symbol_type=symbol_type)
117 return spg
120def get_primitive_structure(
121 structure: Atoms,
122 no_idealize: bool = True,
123 to_primitive: bool = True,
124 symprec: float = 1e-5,
125) -> Atoms:
126 """Returns the primitive structure using spglib.
127 This is a convenience interface to the :func:`standardize_cell`
128 function of spglib that works directly with ase Atoms objects.
130 Parameters
131 ----------
132 structure
133 Input atomic structure.
134 no_idealize
135 If ``True`` lengths and angles are not idealized.
136 to_primitive
137 If ``True`` convert to primitive structure.
138 symprec
139 Tolerance imposed when analyzing the symmetry.
140 """
141 if not spglib_available:
142 raise ImportError('\
143 spglib must be available in order for this function to work, \
144 or you are using a too old version (less than 2.7).') # pragma: no cover
146 structure_tuple = (
147 structure.get_cell(),
148 structure.get_scaled_positions(),
149 structure.numbers)
150 try:
151 result = spglib.standardize_cell(
152 structure_tuple, to_primitive=to_primitive,
153 no_idealize=no_idealize, symprec=symprec)
154 except SpglibCppError:
155 result = None
156 if result is None:
157 raise ValueError('spglib failed to find the primitive cell, maybe caused by large symprec.')
158 lattice, scaled_positions, numbers = result
159 scaled_positions = [np.round(pos, 12) for pos in scaled_positions]
160 structure_prim = Atoms(scaled_positions=scaled_positions,
161 numbers=numbers, cell=lattice, pbc=True)
162 structure_prim.wrap()
164 return structure_prim
167def get_wyckoff_sites(
168 structure: Atoms,
169 map_occupations: List[List[str]] = None,
170 symprec: float = 1e-5,
171 include_representative_atom_index: bool = False,
172) -> List[str]:
173 """Returns the Wyckoff symbols of the input structure.
174 The Wyckoff labels can be conveniently attached as an array to the
175 structure object as demonstrated in the examples section below.
177 By default the occupation of the sites is part of the symmetry
178 analysis. If a chemically disordered structure is provided this
179 will usually reduce the symmetry substantially. If one is
180 interested in the symmetry of the underlying structure one can
181 control how occupations are handled. To this end, one can provide
182 the :attr:`map_occupations` keyword argument. The latter must be a
183 list, each entry of which is a list of species that should be
184 treated as indistinguishable. As a shortcut, if *all* species
185 should be treated as indistinguishable one can provide an empty
186 list. Examples that illustrate the usage of the keyword are given
187 below.
189 Parameters
190 ----------
191 structure
192 Input structure. Note that the occupation of the sites is
193 included in the symmetry analysis.
194 map_occupations
195 Each sublist in this list specifies a group of chemical
196 species that shall be treated as indistinguishable for the
197 purpose of the symmetry analysis.
198 symprec
199 Tolerance imposed when analyzing the symmetry using spglib.
200 include_representative_atom_index
201 If True the index of the first atom in the structure that is
202 representative of the Wyckoff site is included in the symbol.
203 This is in particular useful in cases when there are multiple
204 Wyckoff sites sites with the same Wyckoff letter.
206 Examples
207 --------
208 Wyckoff sites of a hexagonal-close packed structure::
210 >>> from ase.build import bulk
211 >>> structure = bulk('Ti')
212 >>> wyckoff_sites = get_wyckoff_sites(structure)
213 >>> print(wyckoff_sites)
214 ['2d', '2d']
217 The Wyckoff labels can also be attached as an array to the
218 structure, in which case the information is also included when
219 storing the Atoms object::
221 >>> from ase.io import write
222 >>> structure.new_array('wyckoff_sites', wyckoff_sites, str)
223 >>> write('structure.xyz', structure)
225 The function can also be applied to supercells::
227 >>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2)
228 >>> wyckoff_sites = get_wyckoff_sites(structure)
229 >>> print(wyckoff_sites)
230 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
231 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
233 Now assume that one is given a supercell of a (Ga,Al)As
234 alloy. Applying the function directly yields much lower symmetry
235 since the symmetry of the original structure is broken::
237 >>> structure.set_chemical_symbols(
238 ... ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As',
239 ... 'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As'])
240 >>> print(get_wyckoff_sites(structure))
241 ['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i',
242 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i']
244 Since Ga and Al occupy the same sublattice, they should, however,
245 be treated as indistinguishable for the purpose of the symmetry
246 analysis, which can be achieved via the :attr:`map_occupations`
247 keyword::
249 >>> print(get_wyckoff_sites(structure, map_occupations=[['Ga', 'Al'], ['As']]))
250 ['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
251 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']
253 If occupations are to ignored entirely, one can simply provide an
254 empty list. In the present case, this turns the zincblende lattice
255 into a diamond lattice, on which case there is only one Wyckoff
256 site::
258 >>> print(get_wyckoff_sites(structure, map_occupations=[]))
259 ['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a',
260 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a']
261 """
262 if not spglib_available:
263 raise ImportError('\
264 spglib must be available in order for this function to work, \
265 or you are using a too old version (less than 2.7).') # pragma: no cover
267 structure_copy = structure.copy()
268 if map_occupations is not None:
269 if len(map_occupations) > 0:
270 new_symbols = []
271 for symb in structure_copy.get_chemical_symbols():
272 for group in map_occupations: # pragma: no cover - because of the break
273 if symb in group:
274 new_symbols.append(group[0])
275 break
276 else:
277 new_symbols = len(structure) * ['H']
278 structure_copy.set_chemical_symbols(new_symbols)
279 structure_tuple = (
280 structure_copy.get_cell(),
281 structure_copy.get_scaled_positions(),
282 structure_copy.numbers)
283 dataset = spglib.get_symmetry_dataset(structure_tuple, symprec=symprec)
284 n_unitcells = np.linalg.det(dataset.transformation_matrix)
286 equivalent_atoms = list(dataset.equivalent_atoms)
287 wyckoffs = {}
288 for index in set(equivalent_atoms):
289 multiplicity = list(dataset.equivalent_atoms).count(index) / n_unitcells
290 multiplicity = int(round(multiplicity))
291 wyckoff = '{}{}'.format(multiplicity, dataset.wyckoffs[index])
292 if include_representative_atom_index:
293 wyckoff += f'-{index}'
294 wyckoffs[index] = wyckoff
296 return [wyckoffs[equivalent_atoms[a.index]] for a in structure_copy]