Coverage for icet/tools/constituent_strain.py: 99%
147 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1import numpy as np
2import itertools
3from ase import Atoms
4from typing import List, Tuple, Callable
5from ase.symbols import chemical_symbols as ase_chemical_symbols
6from functools import partial
7from .constituent_strain_helper_functions import (_get_structure_factor,
8 _get_partial_structure_factor)
11class KPoint:
12 """
13 Class for handling each k point in a supercell separately.
15 Parameters
16 ----------
17 kpt
18 k-point coordinates.
19 multiplicity
20 Multiplicity of this k-point.
21 structure_factor
22 Current structure associated with this k-point.
23 strain_energy_function
24 Function that takes a concentration and a list of parameters
25 and returns strain energy.
26 damping
27 Damping at this k-point in units of Ångstrom.
28 """
29 def __init__(self, kpt: np.ndarray, multiplicity: float, structure_factor: float,
30 strain_energy_function: Callable[[float, List[float]], float],
31 damping: float):
33 self.kpt = kpt
34 self.multiplicity = multiplicity
35 self.structure_factor = structure_factor
36 self.structure_factor_after = self.structure_factor
37 self.strain_energy_function = strain_energy_function
38 self.damping_factor = np.exp(-(damping * np.linalg.norm(self.kpt)) ** 2)
41class ConstituentStrain:
42 r"""
43 Class for handling constituent strain in cluster expansions
44 (see Laks et al., Phys. Rev. B **46**, 12587 (1992) [LakFerFro92]_).
45 This makes it possible to use cluster expansions to describe systems
46 with strain due to, for example, coherent phase separation.
47 For an extensive example on how to use this module,
48 please see :ref:`this example <constituent_strain_example>`.
50 Parameters
51 ----------
52 supercell
53 Defines supercell that will be used when
54 calculating constituent strain.
55 primitive_structure
56 Primitive structure the supercell is based on.
57 chemical_symbols
58 List with chemical symbols involved, such as ``['Ag', 'Cu']``.
59 concentration_symbol
60 Chemical symbol used to define concentration,
61 such as ``'Ag'``.
62 strain_energy_function
63 A function that takes two arguments, a list of parameters and
64 concentration (e.g., ``[0.5, 0.5, 0.5]`` and ``0.3``),
65 and returns the corresponding strain energy.
66 The parameters are in turn determined by ``k_to_parameter_function``
67 (see below). If ``k_to_parameter_function`` is None,
68 the parameters list will be the k-point. For more information, see
69 :ref:`this example <constituent_strain_example>`.
70 k_to_parameter_function
71 A function that takes a k-point as a list of three floats
72 and returns a parameter vector that will be fed into
73 :attr:`strain_energy_function` (see above). If ``None``, the k-point
74 itself will be the parameter vector to :attr:`strain_energy_function`.
75 The purpose of this function is to be able to precompute
76 any factor in the strain energy that depends on the k-point
77 but not the concentration. For more information, see
78 :ref:`this example <constituent_strain_example>`.
79 damping
80 Damping factor :math:`\eta` used to suppress impact of
81 large-magnitude k-points by multiplying strain with
82 :math:`\exp(-(\eta \mathbf{k})^2)` (unit Ångstrom).
83 tol
84 Numerical tolerance when comparing k-points (units of
85 inverse Ångstrom).
86 """
88 def __init__(self,
89 supercell: Atoms,
90 primitive_structure: Atoms,
91 chemical_symbols: List[str],
92 concentration_symbol: str,
93 strain_energy_function: Callable[[float, List[float]], float],
94 k_to_parameter_function: Callable[[List[float]], List[float]] = None,
95 damping: float = 1.0,
96 tol: float = 1e-6):
97 self.natoms = len(supercell)
99 if len(chemical_symbols) < 2:
100 raise ValueError('Please specify two chemical symbols.')
101 elif len(chemical_symbols) > 2:
102 raise NotImplementedError('The constituent strain module currently'
103 ' only works for binary systems.')
104 spins = [1, -1]
105 self.spin_variables = {ase_chemical_symbols.index(sym): spin
106 for sym, spin in zip(sorted(chemical_symbols), spins)}
107 for key, value in self.spin_variables.items(): 107 ↛ 112line 107 didn't jump to line 112, because the loop on line 107 didn't complete
108 if value == 1: 108 ↛ 107line 108 didn't jump to line 107, because the condition on line 108 was never false
109 self.spin_up = key
110 break
112 self.supercell = supercell
113 self.concentration_number = ase_chemical_symbols.index(concentration_symbol)
115 # Initialize k-points for this supercell
116 self.kpoints = []
117 initial_occupations = self.supercell.get_atomic_numbers()
118 for kpt, multiplicity in _generate_k_points(primitive_structure, supercell, tol=tol):
119 if np.allclose(kpt, [0, 0, 0], atol=tol):
120 continue
122 S = _get_structure_factor(occupations=initial_occupations,
123 positions=self.supercell.get_positions(),
124 kpt=kpt,
125 spin_up=self.spin_up)
127 if k_to_parameter_function is None:
128 parameters = kpt
129 else:
130 parameters = k_to_parameter_function(kpt)
131 self.kpoints.append(KPoint(kpt=kpt,
132 multiplicity=multiplicity,
133 structure_factor=S,
134 damping=damping,
135 strain_energy_function=partial(strain_energy_function,
136 parameters)))
138 def get_concentration(self, occupations: np.ndarray) -> float:
139 """
140 Calculate current concentration.
142 Parameters
143 ----------
144 occupations
145 Current occupations.
146 """
147 return sum(occupations == self.concentration_number) / len(occupations)
149 def _get_constituent_strain_term(self, kpoint: KPoint,
150 occupations: List[int],
151 concentration: float) -> float:
152 """Calculate constituent strain corresponding to a specific k-point.
153 That value is returned and also stored in the corresponding
154 :class:`KPoint <icet.tools.constituent_strain.KPoint>`.
156 Parameters
157 ----------
158 kpoint
159 The k-point to be calculated.
160 occupations
161 Current occupations of the structure.
162 concentration
163 Concentration in the structure.
165 """
166 if abs(concentration) < 1e-9 or abs(1 - concentration) < 1e-9:
167 kpoint.structure_factor = 0.0
168 return 0.0
169 else:
170 # Calculate structure factor
171 S = _get_structure_factor(
172 occupations, self.supercell.get_positions(), kpoint.kpt, self.spin_up)
174 # Save it for faster calculation later on
175 kpoint.structure_factor = S
177 # Constituent strain excluding structure factor
178 DE_CS = kpoint.strain_energy_function(concentration)
179 return DE_CS * kpoint.multiplicity * abs(S)**2 * kpoint.damping_factor / \
180 (4 * concentration * (1 - concentration))
182 def get_constituent_strain(self, occupations: List[int]) -> float:
183 """
184 Calculate total constituent strain.
186 Parameters
187 ----------
188 occupations
189 Current occupations.
190 """
191 c = self.get_concentration(occupations)
192 E_CS = 0.0
193 for kpoint in self.kpoints:
194 E_CS += self._get_constituent_strain_term(kpoint=kpoint,
195 occupations=occupations,
196 concentration=c)
197 return E_CS
199 def get_constituent_strain_change(self, occupations: np.ndarray, atom_index: int) -> float:
200 """
201 Calculate change in constituent strain upon change of the
202 occupation of one site.
204 .. warning ::
205 This function is dependent on the internal state of the
206 :class:`ConstituentStrain` object and **should typically only
207 be used internally by mchammer**. Specifically, the structure
208 factor is saved internally to speed up computation.
209 The first time this function is called, :attr:`occupations`
210 must be the same array as was used to initialize the
211 :class:`ConstituentStrain` object, or the same as was last used
212 when :func:`get_constituent_strain` was called. After the
213 present function has been called, the same occupations vector
214 need to be used the next time as well, unless :func:`accept_change`
215 has been called, in which case :attr:`occupations` should incorporate
216 the changes implied by the previous call to the function.
218 Parameters
219 ----------
220 occupations
221 Occupations before change.
222 atom_index
223 Index of site the occupation of which is to be changed.
224 """
225 # Determine concentration before and after
226 n = sum(occupations == self.concentration_number)
227 c_before = n / self.natoms
228 if occupations[atom_index] == self.concentration_number:
229 c_after = (n - 1) / self.natoms
230 else:
231 c_after = (n + 1) / self.natoms
233 E_CS_change = 0.0
234 position = self.supercell[atom_index].position
235 spin = self.spin_variables[occupations[atom_index]]
236 for kpoint in self.kpoints:
237 # Change in structure factor
238 S_before = kpoint.structure_factor
239 dS = _get_partial_structure_factor(kpoint.kpt, position, self.natoms)
240 S_after = S_before - 2 * spin * dS
242 # Save the new value in a variable such that it can be
243 # accessed later if the change is accepted
244 kpoint.structure_factor_after = S_after
246 # Energy before
247 if abs(c_before) < 1e-9 or abs(c_before - 1) < 1e-9:
248 E_before = 0.0
249 else:
250 E_before = kpoint.strain_energy_function(c_before)
251 E_before *= abs(S_before)**2 / (4 * c_before * (1 - c_before))
253 # Energy after
254 if abs(c_after) < 1e-9 or abs(c_after - 1) < 1e-9:
255 E_after = 0.0
256 else:
257 E_after = kpoint.strain_energy_function(c_after)
258 E_after *= abs(S_after)**2 / (4 * c_after * (1 - c_after))
260 # Difference
261 E_CS_change += kpoint.multiplicity * kpoint.damping_factor * (E_after - E_before)
262 return E_CS_change
264 def accept_change(self) -> None:
265 """
266 Update structure factor for each kpoint to the value
267 in :attr:`structure_factor_after`. This makes it possible to
268 efficiently calculate changes in constituent strain with
269 the :func:`get_constituent_strain_change` function; this function
270 should be called if the last occupations used to call
271 :func:`get_constituent_strain_change` should be the starting point
272 for the next call of :func:`get_constituent_strain_change`.
273 This is taken care of automatically by the Monte Carlo
274 simulations in :program:`mchammer`.
275 """
276 for kpoint in self.kpoints:
277 kpoint.structure_factor = kpoint.structure_factor_after
280def _generate_k_points(primitive_structure: Atoms,
281 supercell: Atoms,
282 tol: float) -> Tuple[np.ndarray, float]:
283 """
284 Generate all k-points in the 1BZ of the primitive cell
285 that potentially correspond to a nonzero structure factor.
286 These are all the k-points that are integer multiples of
287 the reciprocal supercell.
289 Parameters
290 ----------
291 primitive_structure
292 Primitive structure that the supercell is based on.
293 supercell
294 Supercell of primitive structure.
295 """
296 reciprocal_primitive = np.linalg.inv(primitive_structure.cell) # column vectors
297 reciprocal_supercell = np.linalg.inv(supercell.cell) # column vectors
299 # How many k-points are there?
300 # This information is needed to know when to stop looking for more,
301 # i.e., it implicitly determines the iteration limits.
302 nkpoints = int(round(np.linalg.det(reciprocal_primitive)
303 / np.linalg.det(reciprocal_supercell)))
305 # The gamma point is always present
306 yield np.array([0., 0., 0.]), 1.0
307 found_kpoints = [np.array([0., 0., 0.])]
308 covered_nkpoints = 1
310 # Now loop until we have found all k-points.
311 # We loop by successively increasing the sum of the
312 # absolute values of the components of the reciprocal
313 # supercell vectors, and we stop when we have
314 # found the correct number of k-points.
315 component_sum = 0
316 while covered_nkpoints < nkpoints: 316 ↛ exitline 316 didn't return from function '_generate_k_points', because the condition on line 316 was never false
317 component_sum += 1
318 for ijk_p in _ordered_combinations(component_sum, 3):
319 for signs in itertools.product([-1, 1], repeat=3):
320 if np.any([ijk_p[i] == 0 and signs[i] < 0 for i in range(3)]):
321 continue
322 q_sc = np.multiply(ijk_p, signs)
323 k = np.dot(reciprocal_supercell, q_sc)
325 # Now check whether we can come closer to the gamma point
326 # by translation with primitive reciprocal lattice vectors.
327 # In other words, we translate the point to the 1BZ.
328 k = _translate_to_1BZ(k, reciprocal_primitive, tol=tol)
329 equivalent_kpoints = _find_equivalent_kpoints(k, reciprocal_primitive, tol=tol)
331 # Have we already found this k-point?
332 found = False
333 for k in equivalent_kpoints:
334 for k_comp in found_kpoints:
335 if np.linalg.norm(k - k_comp) < tol:
336 # Then we had already found it
337 found = True
338 break
339 if found:
340 break
341 else:
342 # Then this is a new k-point
343 # Yield it and its equivalent friends
344 covered_nkpoints += 1
345 for k in equivalent_kpoints:
346 found_kpoints.append(k)
347 yield k, 1 / len(equivalent_kpoints)
349 # Break if we have found all k-points
350 assert covered_nkpoints <= nkpoints
351 if covered_nkpoints == nkpoints:
352 break
355def _ordered_combinations(s: int, n: int) -> List[int]:
356 """
357 Recursively generate all combinations of n integers
358 that sum to s.
360 Parameters
361 ----------
362 s
363 Required sum of the combination
364 n
365 Number of intergers in the list
366 """
367 if n == 1:
368 yield [s]
369 else:
370 for i in range(s + 1):
371 for rest in _ordered_combinations(s - i, n - 1):
372 rest.append(i)
373 yield rest
376def _translate_to_1BZ(kpt: np.ndarray, primitive: np.ndarray, tol: float) -> np.ndarray:
377 """
378 Translate k-point into 1BZ by translating it by
379 primitive lattice vectors and testing if that takes
380 us closer to gamma. The algorithm works by recursion
381 such that we keep translating until we no longer get
382 closer (will this always work?).
384 Parameters
385 ----------
386 kpt
387 coordinates of k-point that we translate
388 primitive
389 primitive reciprocal cell with lattice vectors as
390 column vectors
391 tol
392 numerical tolerance when comparing k-points (units of
393 inverse Ångstrom)
395 Returns
396 -------
397 the k-point translated into 1BZ
398 """
399 original_distance_from_gamma = np.linalg.norm(kpt)
400 min_distance_from_gamma = original_distance_from_gamma
401 best_kpt = kpt
402 # Loop through translations and see if any takes us closer to gamma
403 for translation in _generate_primitive_translations(primitive):
404 new_kpt = kpt + translation
405 new_distance_from_gamma = np.linalg.norm(new_kpt)
406 if new_distance_from_gamma < min_distance_from_gamma:
407 min_distance_from_gamma = new_distance_from_gamma
408 best_kpt = new_kpt
409 if abs(original_distance_from_gamma - min_distance_from_gamma) < tol:
410 # Then we did not come closer to gamma, so we are done
411 return best_kpt
412 else:
413 # We came closer to gamma. Maybe we can come even closer?
414 return _translate_to_1BZ(best_kpt, primitive, tol=tol)
417def _find_equivalent_kpoints(kpt: np.ndarray,
418 primitive: np.ndarray,
419 tol: float) -> List[np.ndarray]:
420 """
421 Find k-points in 1BZ equivalent to this k-point, i.e.,
422 points in the 1BZ related by a primitive translation
423 (for example, if a k-point lies at the edge of the 1BZ, we
424 will find it on the opposite side of the 1BZ too).
426 Parameters
427 ----------
428 kpt
429 k-point coordinates
430 primitive
431 primitive reciprocal lattice vectors as columns
432 tol
433 numerical tolerance when comparing k-points (units of
434 inverse Ångstrom)
436 Returns
437 -------
438 list of equivalent k-points
439 """
440 original_norm = np.linalg.norm(kpt)
441 equivalent_kpoints = [kpt]
442 for translation in _generate_primitive_translations(primitive):
443 if abs(original_norm - np.linalg.norm(kpt + translation)) < tol:
444 equivalent_kpoints.append(kpt + translation)
445 return equivalent_kpoints
448def _generate_primitive_translations(primitive: np.ndarray) -> np.ndarray:
449 """
450 Generate the 26 (3x3x3 - 1) primitive translation vectors
451 defined by (0/+1/-1, 0/+1/-1, 0/+1/-1) (exlcluding (0, 0, 0))
453 Parameters
454 ----------
455 primitive
456 Lattice vectors as columns
458 Yields
459 ------
460 traslation vector, a multiple of primitive lattice vectors
461 """
462 for ijk in itertools.product((0, 1, -1), repeat=3):
463 if ijk == (0, 0, 0):
464 continue
465 yield np.dot(primitive, ijk)