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