Coverage for /builds/ase/ase/ase/ga/standardmutations.py: 74.06%
374 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""A collection of mutations that can be used."""
4from math import cos, pi, sin
6import numpy as np
8from ase import Atoms
9from ase.calculators.lammps.coordinatetransform import calc_rotated_cell
10from ase.cell import Cell
11from ase.ga.offspring_creator import CombinationMutation, OffspringCreator
12from ase.ga.utilities import (
13 atoms_too_close,
14 atoms_too_close_two_sets,
15 gather_atoms_by_tag,
16 get_rotation_matrix,
17)
20class RattleMutation(OffspringCreator):
21 """An implementation of the rattle mutation as described in:
23 R.L. Johnston Dalton Transactions, Vol. 22,
24 No. 22. (2003), pp. 4193-4207
26 Parameters:
28 blmin: Dictionary defining the minimum distance between atoms
29 after the rattle.
31 n_top: Number of atoms optimized by the GA.
33 rattle_strength: Strength with which the atoms are moved.
35 rattle_prop: The probability with which each atom is rattled.
37 test_dist_to_slab: whether to also make sure that the distances
38 between the atoms and the slab satisfy the blmin.
40 use_tags: if True, the atomic tags will be used to preserve
41 molecular identity. Same-tag atoms will then be
42 displaced collectively, so that the internal
43 geometry is preserved.
45 rng: Random number generator
46 By default numpy.random.
47 """
49 def __init__(self, blmin, n_top, rattle_strength=0.8,
50 rattle_prop=0.4, test_dist_to_slab=True, use_tags=False,
51 verbose=False, rng=np.random):
52 OffspringCreator.__init__(self, verbose, rng=rng)
53 self.blmin = blmin
54 self.n_top = n_top
55 self.rattle_strength = rattle_strength
56 self.rattle_prop = rattle_prop
57 self.test_dist_to_slab = test_dist_to_slab
58 self.use_tags = use_tags
60 self.descriptor = 'RattleMutation'
61 self.min_inputs = 1
63 def get_new_individual(self, parents):
64 f = parents[0]
66 indi = self.mutate(f)
67 if indi is None:
68 return indi, 'mutation: rattle'
70 indi = self.initialize_individual(f, indi)
71 indi.info['data']['parents'] = [f.info['confid']]
73 return self.finalize_individual(indi), 'mutation: rattle'
75 def mutate(self, atoms):
76 """Does the actual mutation."""
77 N = len(atoms) if self.n_top is None else self.n_top
78 slab = atoms[:len(atoms) - N]
79 atoms = atoms[-N:]
80 tags = atoms.get_tags() if self.use_tags else np.arange(N)
81 pos_ref = atoms.get_positions()
82 num = atoms.get_atomic_numbers()
83 cell = atoms.get_cell()
84 pbc = atoms.get_pbc()
85 st = 2. * self.rattle_strength
87 count = 0
88 maxcount = 1000
89 too_close = True
90 while too_close and count < maxcount:
91 count += 1
92 pos = pos_ref.copy()
93 ok = False
94 for tag in np.unique(tags):
95 select = np.where(tags == tag)
96 if self.rng.random() < self.rattle_prop:
97 ok = True
98 r = self.rng.random(3)
99 pos[select] += st * (r - 0.5)
101 if not ok:
102 # Nothing got rattled
103 continue
105 top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
106 too_close = atoms_too_close(
107 top, self.blmin, use_tags=self.use_tags)
108 if not too_close and self.test_dist_to_slab:
109 too_close = atoms_too_close_two_sets(top, slab, self.blmin)
111 if count == maxcount:
112 return None
114 mutant = slab + top
115 return mutant
118class PermutationMutation(OffspringCreator):
119 """Mutation that permutes a percentage of the atom types in the cluster.
121 Parameters:
123 n_top: Number of atoms optimized by the GA.
125 probability: The probability with which an atom is permuted.
127 test_dist_to_slab: whether to also make sure that the distances
128 between the atoms and the slab satisfy the blmin.
130 use_tags: if True, the atomic tags will be used to preserve
131 molecular identity. Permutations will then happen
132 at the molecular level, i.e. swapping the center-of-
133 positions of two moieties while preserving their
134 internal geometries.
136 blmin: Dictionary defining the minimum distance between atoms
137 after the permutation. If equal to None (the default),
138 no such check is performed.
140 rng: Random number generator
141 By default numpy.random.
142 """
144 def __init__(self, n_top, probability=0.33, test_dist_to_slab=True,
145 use_tags=False, blmin=None, rng=np.random, verbose=False):
146 OffspringCreator.__init__(self, verbose, rng=rng)
147 self.n_top = n_top
148 self.probability = probability
149 self.test_dist_to_slab = test_dist_to_slab
150 self.use_tags = use_tags
151 self.blmin = blmin
153 self.descriptor = 'PermutationMutation'
154 self.min_inputs = 1
156 def get_new_individual(self, parents):
157 f = parents[0]
159 indi = self.mutate(f)
160 if indi is None:
161 return indi, 'mutation: permutation'
163 indi = self.initialize_individual(f, indi)
164 indi.info['data']['parents'] = [f.info['confid']]
166 return self.finalize_individual(indi), 'mutation: permutation'
168 def mutate(self, atoms):
169 """Does the actual mutation."""
170 N = len(atoms) if self.n_top is None else self.n_top
171 slab = atoms[:len(atoms) - N]
172 atoms = atoms[-N:]
173 if self.use_tags:
174 gather_atoms_by_tag(atoms)
175 tags = atoms.get_tags() if self.use_tags else np.arange(N)
176 pos_ref = atoms.get_positions()
177 num = atoms.get_atomic_numbers()
178 cell = atoms.get_cell()
179 pbc = atoms.get_pbc()
180 symbols = atoms.get_chemical_symbols()
182 unique_tags = np.unique(tags)
183 n = len(unique_tags)
184 swaps = int(np.ceil(n * self.probability / 2.))
186 sym = []
187 for tag in unique_tags:
188 indices = np.where(tags == tag)[0]
189 s = ''.join([symbols[j] for j in indices])
190 sym.append(s)
192 assert len(np.unique(sym)) > 1, \
193 'Permutations with one atom (or molecule) type is not valid'
195 count = 0
196 maxcount = 1000
197 too_close = True
198 while too_close and count < maxcount:
199 count += 1
200 pos = pos_ref.copy()
201 for _ in range(swaps):
202 i = j = 0
203 while sym[i] == sym[j]:
204 i = self.rng.randint(0, high=n)
205 j = self.rng.randint(0, high=n)
206 ind1 = np.where(tags == i)
207 ind2 = np.where(tags == j)
208 cop1 = np.mean(pos[ind1], axis=0)
209 cop2 = np.mean(pos[ind2], axis=0)
210 pos[ind1] += cop2 - cop1
211 pos[ind2] += cop1 - cop2
213 top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
214 if self.blmin is None:
215 too_close = False
216 else:
217 too_close = atoms_too_close(
218 top, self.blmin, use_tags=self.use_tags)
219 if not too_close and self.test_dist_to_slab:
220 too_close = atoms_too_close_two_sets(top, slab, self.blmin)
222 if count == maxcount:
223 return None
225 mutant = slab + top
226 return mutant
229class MirrorMutation(OffspringCreator):
230 """A mirror mutation, as described in
231 TO BE PUBLISHED.
233 This mutation mirrors half of the cluster in a
234 randomly oriented cutting plane discarding the other half.
236 Parameters:
238 blmin: Dictionary defining the minimum allowed
239 distance between atoms.
241 n_top: Number of atoms the GA optimizes.
243 reflect: Defines if the mirrored half is also reflected
244 perpendicular to the mirroring plane.
246 rng: Random number generator
247 By default numpy.random.
248 """
250 def __init__(self, blmin, n_top, reflect=False, rng=np.random,
251 verbose=False):
252 OffspringCreator.__init__(self, verbose, rng=rng)
253 self.blmin = blmin
254 self.n_top = n_top
255 self.reflect = reflect
257 self.descriptor = 'MirrorMutation'
258 self.min_inputs = 1
260 def get_new_individual(self, parents):
261 f = parents[0]
263 indi = self.mutate(f)
264 if indi is None:
265 return indi, 'mutation: mirror'
267 indi = self.initialize_individual(f, indi)
268 indi.info['data']['parents'] = [f.info['confid']]
270 return self.finalize_individual(indi), 'mutation: mirror'
272 def mutate(self, atoms):
273 """ Do the mutation of the atoms input. """
274 reflect = self.reflect
275 tc = True
276 slab = atoms[0:len(atoms) - self.n_top]
277 top = atoms[len(atoms) - self.n_top: len(atoms)]
278 num = top.numbers
279 unique_types = list(set(num))
280 nu = {u: sum(num == u) for u in unique_types}
281 n_tries = 1000
282 counter = 0
283 changed = False
285 while tc and counter < n_tries:
286 counter += 1
287 cand = top.copy()
288 pos = cand.get_positions()
290 cm = np.average(top.get_positions(), axis=0)
292 # first select a randomly oriented cutting plane
293 theta = pi * self.rng.random()
294 phi = 2. * pi * self.rng.random()
295 n = (cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta))
296 n = np.array(n)
298 # Calculate all atoms signed distance to the cutting plane
299 D = []
300 for (i, p) in enumerate(pos):
301 d = np.dot(p - cm, n)
302 D.append((i, d))
304 # Sort the atoms by their signed distance
305 D.sort(key=lambda x: x[1])
306 nu_taken = {}
308 # Select half of the atoms needed for a full cluster
309 p_use = []
310 n_use = []
311 for (i, d) in D:
312 if num[i] not in nu_taken.keys():
313 nu_taken[num[i]] = 0
314 if nu_taken[num[i]] < nu[num[i]] / 2.:
315 p_use.append(pos[i])
316 n_use.append(num[i])
317 nu_taken[num[i]] += 1
319 # calculate the mirrored position and add these.
320 pn = []
321 for p in p_use:
322 pt = p - 2. * np.dot(p - cm, n) * n
323 if reflect:
324 pt = -pt + 2 * cm + 2 * n * np.dot(pt - cm, n)
325 pn.append(pt)
327 n_use.extend(n_use)
328 p_use.extend(pn)
330 # In the case of an uneven number of
331 # atoms we need to add one extra
332 for n in nu:
333 if nu[n] % 2 == 0:
334 continue
335 while sum(n_use == n) > nu[n]:
336 for i in range(int(len(n_use) / 2), len(n_use)):
337 if n_use[i] == n:
338 del p_use[i]
339 del n_use[i]
340 break
341 assert sum(n_use == n) == nu[n]
343 # Make sure we have the correct number of atoms
344 # and rearrange the atoms so they are in the right order
345 for i in range(len(n_use)):
346 if num[i] == n_use[i]:
347 continue
348 for j in range(i + 1, len(n_use)):
349 if n_use[j] == num[i]:
350 tn = n_use[i]
351 tp = p_use[i]
352 n_use[i] = n_use[j]
353 p_use[i] = p_use[j]
354 p_use[j] = tp
355 n_use[j] = tn
357 # Finally we check that nothing is too close in the end product.
358 cand = Atoms(num, p_use, cell=slab.get_cell(), pbc=slab.get_pbc())
360 tc = atoms_too_close(cand, self.blmin)
361 if tc:
362 continue
363 tc = atoms_too_close_two_sets(slab, cand, self.blmin)
365 if not changed and counter > n_tries // 2:
366 reflect = not reflect
367 changed = True
369 tot = slab + cand
371 if counter == n_tries:
372 return None
373 return tot
376class StrainMutation(OffspringCreator):
377 """ Mutates a candidate by applying a randomly generated strain.
379 For more information, see also:
381 * :doi:`Glass, Oganov, Hansen, Comp. Phys. Comm. 175 (2006) 713-720
382 <10.1016/j.cpc.2006.07.020>`
384 * :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
385 <10.1016/j.cpc.2010.07.048>`
387 After initialization of the mutation, a scaling volume
388 (to which each mutated structure is scaled before checking the
389 constraints) is typically generated from the population,
390 which is then also occasionally updated in the course of the
391 GA run.
393 Parameters:
395 blmin: dict
396 The closest allowed interatomic distances on the form:
397 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
399 cellbounds: ase.ga.utilities.CellBounds instance
400 Describes limits on the cell shape, see
401 :class:`~ase.ga.utilities.CellBounds`.
403 stddev: float
404 Standard deviation used in the generation of the
405 strain matrix elements.
407 number_of_variable_cell_vectors: int (default 3)
408 The number of variable cell vectors (1, 2 or 3).
409 To keep things simple, it is the 'first' vectors which
410 will be treated as variable, i.e. the 'a' vector in the
411 univariate case, the 'a' and 'b' vectors in the bivariate
412 case, etc.
414 use_tags: boolean
415 Whether to use the atomic tags to preserve molecular identity.
417 rng: Random number generator
418 By default numpy.random.
419 """
421 def __init__(self, blmin, cellbounds=None, stddev=0.7,
422 number_of_variable_cell_vectors=3, use_tags=False,
423 rng=np.random, verbose=False):
424 OffspringCreator.__init__(self, verbose, rng=rng)
425 self.blmin = blmin
426 self.cellbounds = cellbounds
427 self.stddev = stddev
428 self.number_of_variable_cell_vectors = number_of_variable_cell_vectors
429 self.use_tags = use_tags
431 self.scaling_volume = None
432 self.descriptor = 'StrainMutation'
433 self.min_inputs = 1
435 def update_scaling_volume(self, population, w_adapt=0.5, n_adapt=0):
436 """Function to initialize or update the scaling volume in a GA run.
438 w_adapt: weight of the new vs the old scaling volume
440 n_adapt: number of best candidates in the population that
441 are used to calculate the new scaling volume
442 """
443 if not n_adapt:
444 # if not set, take best 20% of the population
445 n_adapt = int(np.ceil(0.2 * len(population)))
446 v_new = np.mean([a.get_volume() for a in population[:n_adapt]])
448 if not self.scaling_volume:
449 self.scaling_volume = v_new
450 else:
451 volumes = [self.scaling_volume, v_new]
452 weights = [1 - w_adapt, w_adapt]
453 self.scaling_volume = np.average(volumes, weights=weights)
455 def get_new_individual(self, parents):
456 f = parents[0]
458 indi = self.mutate(f)
459 if indi is None:
460 return indi, 'mutation: strain'
462 indi = self.initialize_individual(f, indi)
463 indi.info['data']['parents'] = [f.info['confid']]
465 return self.finalize_individual(indi), 'mutation: strain'
467 def mutate(self, atoms):
468 """ Does the actual mutation. """
469 cell_ref = atoms.get_cell()
470 pos_ref = atoms.get_positions()
472 if self.scaling_volume is None:
473 # The scaling_volume has not been set (yet),
474 # so we give it the same volume as the parent
475 vol_ref = atoms.get_volume()
476 else:
477 vol_ref = self.scaling_volume
479 if self.use_tags:
480 tags = atoms.get_tags()
481 gather_atoms_by_tag(atoms)
482 pos = atoms.get_positions()
484 mutant = atoms.copy()
486 count = 0
487 too_close = True
488 maxcount = 1000
489 while too_close and count < maxcount:
490 count += 1
492 # generating the strain matrix:
493 strain = np.identity(3)
494 for i in range(self.number_of_variable_cell_vectors):
495 for j in range(i + 1):
496 r = self.rng.normal(loc=0., scale=self.stddev)
497 if i == j:
498 strain[i, j] += r
499 else:
500 epsilon = 0.5 * r
501 strain[i, j] += epsilon
502 strain[j, i] += epsilon
504 # applying the strain:
505 cell_new = np.dot(strain, cell_ref)
507 # convert the submatrix with the variable cell vectors
508 # to a lower triangular form
509 cell_new = calc_rotated_cell(cell_new)
510 for i in range(self.number_of_variable_cell_vectors, 3):
511 cell_new[i] = cell_ref[i]
513 cell_new = Cell(cell_new)
515 # volume scaling:
516 if self.number_of_variable_cell_vectors > 0:
517 scaling = vol_ref / cell_new.volume
518 scaling **= 1. / self.number_of_variable_cell_vectors
519 cell_new[:self.number_of_variable_cell_vectors] *= scaling
521 # check cell dimensions:
522 if not self.cellbounds.is_within_bounds(cell_new):
523 continue
525 # ensure non-variable cell vectors are indeed unchanged
526 for i in range(self.number_of_variable_cell_vectors, 3):
527 assert np.allclose(cell_new[i], cell_ref[i])
529 # check that the volume is correct
530 assert np.isclose(vol_ref, cell_new.volume)
532 # apply the new unit cell and scale
533 # the atomic positions accordingly
534 mutant.set_cell(cell_ref, scale_atoms=False)
536 if self.use_tags:
537 transfo = np.linalg.solve(cell_ref, cell_new)
538 for tag in np.unique(tags):
539 select = np.where(tags == tag)
540 cop = np.mean(pos[select], axis=0)
541 disp = np.dot(cop, transfo) - cop
542 mutant.positions[select] += disp
543 else:
544 mutant.set_positions(pos_ref)
546 mutant.set_cell(cell_new, scale_atoms=not self.use_tags)
547 mutant.wrap()
549 # check the interatomic distances
550 too_close = atoms_too_close(mutant, self.blmin,
551 use_tags=self.use_tags)
553 if count == maxcount:
554 mutant = None
556 return mutant
559class PermuStrainMutation(CombinationMutation):
560 """Combination of PermutationMutation and StrainMutation.
562 For more information, see also:
564 * :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
565 <10.1016/j.cpc.2010.07.048>`
567 Parameters:
569 permutationmutation: OffspringCreator instance
570 A mutation that permutes atom types.
572 strainmutation: OffspringCreator instance
573 A mutation that mutates by straining.
574 """
576 def __init__(self, permutationmutation, strainmutation, verbose=False):
577 super().__init__(permutationmutation,
578 strainmutation,
579 verbose=verbose)
580 self.descriptor = 'permustrain'
583class RotationalMutation(OffspringCreator):
584 """Mutates a candidate by applying random rotations
585 to multi-atom moieties in the structure (atoms with
586 the same tag are considered part of one such moiety).
588 Only performs whole-molecule rotations, no internal
589 rotations.
591 For more information, see also:
593 * `Zhu Q., Oganov A.R., Glass C.W., Stokes H.T,
594 Acta Cryst. (2012), B68, 215-226.`__
596 __ https://dx.doi.org/10.1107/S0108768112017466
598 Parameters:
600 blmin: dict
601 The closest allowed interatomic distances on the form:
602 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
604 n_top: int or None
605 The number of atoms to optimize (None = include all).
607 fraction: float
608 Fraction of the moieties to be rotated.
610 tags: None or list of integers
611 Specifies, respectively, whether all moieties or only those
612 with matching tags are eligible for rotation.
614 min_angle: float
615 Minimal angle (in radians) for each rotation;
616 should lie in the interval [0, pi].
618 test_dist_to_slab: boolean
619 Whether also the distances to the slab
620 should be checked to satisfy the blmin.
622 rng: Random number generator
623 By default numpy.random.
624 """
626 def __init__(self, blmin, n_top=None, fraction=0.33, tags=None,
627 min_angle=1.57, test_dist_to_slab=True, rng=np.random,
628 verbose=False):
629 OffspringCreator.__init__(self, verbose, rng=rng)
630 self.blmin = blmin
631 self.n_top = n_top
632 self.fraction = fraction
633 self.tags = tags
634 self.min_angle = min_angle
635 self.test_dist_to_slab = test_dist_to_slab
636 self.descriptor = 'RotationalMutation'
637 self.min_inputs = 1
639 def get_new_individual(self, parents):
640 f = parents[0]
642 indi = self.mutate(f)
643 if indi is None:
644 return indi, 'mutation: rotational'
646 indi = self.initialize_individual(f, indi)
647 indi.info['data']['parents'] = [f.info['confid']]
649 return self.finalize_individual(indi), 'mutation: rotational'
651 def mutate(self, atoms):
652 """Does the actual mutation."""
653 N = len(atoms) if self.n_top is None else self.n_top
654 slab = atoms[:len(atoms) - N]
655 atoms = atoms[-N:]
657 mutant = atoms.copy()
658 gather_atoms_by_tag(mutant)
659 pos = mutant.get_positions()
660 tags = mutant.get_tags()
661 eligible_tags = tags if self.tags is None else self.tags
663 indices = {}
664 for tag in np.unique(tags):
665 hits = np.where(tags == tag)[0]
666 if len(hits) > 1 and tag in eligible_tags:
667 indices[tag] = hits
669 n_rot = int(np.ceil(len(indices) * self.fraction))
670 chosen_tags = self.rng.choice(list(indices.keys()), size=n_rot,
671 replace=False)
673 too_close = True
674 count = 0
675 maxcount = 10000
676 while too_close and count < maxcount:
677 newpos = np.copy(pos)
678 for tag in chosen_tags:
679 p = np.copy(newpos[indices[tag]])
680 cop = np.mean(p, axis=0)
682 if len(p) == 2:
683 line = (p[1] - p[0]) / np.linalg.norm(p[1] - p[0])
684 while True:
685 axis = self.rng.random(3)
686 axis /= np.linalg.norm(axis)
687 a = np.arccos(np.dot(axis, line))
688 if np.pi / 4 < a < np.pi * 3 / 4:
689 break
690 else:
691 axis = self.rng.random(3)
692 axis /= np.linalg.norm(axis)
694 angle = self.min_angle
695 angle += 2 * (np.pi - self.min_angle) * self.rng.random()
697 m = get_rotation_matrix(axis, angle)
698 newpos[indices[tag]] = np.dot(m, (p - cop).T).T + cop
700 mutant.set_positions(newpos)
701 mutant.wrap()
702 too_close = atoms_too_close(mutant, self.blmin, use_tags=True)
703 count += 1
705 if not too_close and self.test_dist_to_slab:
706 too_close = atoms_too_close_two_sets(slab, mutant, self.blmin)
708 if count == maxcount:
709 mutant = None
710 else:
711 mutant = slab + mutant
713 return mutant
716class RattleRotationalMutation(CombinationMutation):
717 """Combination of RattleMutation and RotationalMutation.
719 Parameters:
721 rattlemutation: OffspringCreator instance
722 A mutation that rattles atoms.
724 rotationalmutation: OffspringCreator instance
725 A mutation that rotates moieties.
726 """
728 def __init__(self, rattlemutation, rotationalmutation, verbose=False):
729 super().__init__(rattlemutation,
730 rotationalmutation,
731 verbose=verbose)
732 self.descriptor = 'rattlerotational'