Coverage for /builds/debichem-team/python-ase/ase/ga/standardmutations.py: 74.06%
374 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1"""A collection of mutations that can be used."""
2from math import cos, pi, sin
4import numpy as np
6from ase import Atoms
7from ase.calculators.lammps.coordinatetransform import calc_rotated_cell
8from ase.cell import Cell
9from ase.ga.offspring_creator import CombinationMutation, OffspringCreator
10from ase.ga.utilities import (
11 atoms_too_close,
12 atoms_too_close_two_sets,
13 gather_atoms_by_tag,
14 get_rotation_matrix,
15)
18class RattleMutation(OffspringCreator):
19 """An implementation of the rattle mutation as described in:
21 R.L. Johnston Dalton Transactions, Vol. 22,
22 No. 22. (2003), pp. 4193-4207
24 Parameters:
26 blmin: Dictionary defining the minimum distance between atoms
27 after the rattle.
29 n_top: Number of atoms optimized by the GA.
31 rattle_strength: Strength with which the atoms are moved.
33 rattle_prop: The probability with which each atom is rattled.
35 test_dist_to_slab: whether to also make sure that the distances
36 between the atoms and the slab satisfy the blmin.
38 use_tags: if True, the atomic tags will be used to preserve
39 molecular identity. Same-tag atoms will then be
40 displaced collectively, so that the internal
41 geometry is preserved.
43 rng: Random number generator
44 By default numpy.random.
45 """
47 def __init__(self, blmin, n_top, rattle_strength=0.8,
48 rattle_prop=0.4, test_dist_to_slab=True, use_tags=False,
49 verbose=False, rng=np.random):
50 OffspringCreator.__init__(self, verbose, rng=rng)
51 self.blmin = blmin
52 self.n_top = n_top
53 self.rattle_strength = rattle_strength
54 self.rattle_prop = rattle_prop
55 self.test_dist_to_slab = test_dist_to_slab
56 self.use_tags = use_tags
58 self.descriptor = 'RattleMutation'
59 self.min_inputs = 1
61 def get_new_individual(self, parents):
62 f = parents[0]
64 indi = self.mutate(f)
65 if indi is None:
66 return indi, 'mutation: rattle'
68 indi = self.initialize_individual(f, indi)
69 indi.info['data']['parents'] = [f.info['confid']]
71 return self.finalize_individual(indi), 'mutation: rattle'
73 def mutate(self, atoms):
74 """Does the actual mutation."""
75 N = len(atoms) if self.n_top is None else self.n_top
76 slab = atoms[:len(atoms) - N]
77 atoms = atoms[-N:]
78 tags = atoms.get_tags() if self.use_tags else np.arange(N)
79 pos_ref = atoms.get_positions()
80 num = atoms.get_atomic_numbers()
81 cell = atoms.get_cell()
82 pbc = atoms.get_pbc()
83 st = 2. * self.rattle_strength
85 count = 0
86 maxcount = 1000
87 too_close = True
88 while too_close and count < maxcount:
89 count += 1
90 pos = pos_ref.copy()
91 ok = False
92 for tag in np.unique(tags):
93 select = np.where(tags == tag)
94 if self.rng.random() < self.rattle_prop:
95 ok = True
96 r = self.rng.random(3)
97 pos[select] += st * (r - 0.5)
99 if not ok:
100 # Nothing got rattled
101 continue
103 top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
104 too_close = atoms_too_close(
105 top, self.blmin, use_tags=self.use_tags)
106 if not too_close and self.test_dist_to_slab:
107 too_close = atoms_too_close_two_sets(top, slab, self.blmin)
109 if count == maxcount:
110 return None
112 mutant = slab + top
113 return mutant
116class PermutationMutation(OffspringCreator):
117 """Mutation that permutes a percentage of the atom types in the cluster.
119 Parameters:
121 n_top: Number of atoms optimized by the GA.
123 probability: The probability with which an atom is permuted.
125 test_dist_to_slab: whether to also make sure that the distances
126 between the atoms and the slab satisfy the blmin.
128 use_tags: if True, the atomic tags will be used to preserve
129 molecular identity. Permutations will then happen
130 at the molecular level, i.e. swapping the center-of-
131 positions of two moieties while preserving their
132 internal geometries.
134 blmin: Dictionary defining the minimum distance between atoms
135 after the permutation. If equal to None (the default),
136 no such check is performed.
138 rng: Random number generator
139 By default numpy.random.
140 """
142 def __init__(self, n_top, probability=0.33, test_dist_to_slab=True,
143 use_tags=False, blmin=None, rng=np.random, verbose=False):
144 OffspringCreator.__init__(self, verbose, rng=rng)
145 self.n_top = n_top
146 self.probability = probability
147 self.test_dist_to_slab = test_dist_to_slab
148 self.use_tags = use_tags
149 self.blmin = blmin
151 self.descriptor = 'PermutationMutation'
152 self.min_inputs = 1
154 def get_new_individual(self, parents):
155 f = parents[0]
157 indi = self.mutate(f)
158 if indi is None:
159 return indi, 'mutation: permutation'
161 indi = self.initialize_individual(f, indi)
162 indi.info['data']['parents'] = [f.info['confid']]
164 return self.finalize_individual(indi), 'mutation: permutation'
166 def mutate(self, atoms):
167 """Does the actual mutation."""
168 N = len(atoms) if self.n_top is None else self.n_top
169 slab = atoms[:len(atoms) - N]
170 atoms = atoms[-N:]
171 if self.use_tags:
172 gather_atoms_by_tag(atoms)
173 tags = atoms.get_tags() if self.use_tags else np.arange(N)
174 pos_ref = atoms.get_positions()
175 num = atoms.get_atomic_numbers()
176 cell = atoms.get_cell()
177 pbc = atoms.get_pbc()
178 symbols = atoms.get_chemical_symbols()
180 unique_tags = np.unique(tags)
181 n = len(unique_tags)
182 swaps = int(np.ceil(n * self.probability / 2.))
184 sym = []
185 for tag in unique_tags:
186 indices = np.where(tags == tag)[0]
187 s = ''.join([symbols[j] for j in indices])
188 sym.append(s)
190 assert len(np.unique(sym)) > 1, \
191 'Permutations with one atom (or molecule) type is not valid'
193 count = 0
194 maxcount = 1000
195 too_close = True
196 while too_close and count < maxcount:
197 count += 1
198 pos = pos_ref.copy()
199 for _ in range(swaps):
200 i = j = 0
201 while sym[i] == sym[j]:
202 i = self.rng.randint(0, high=n)
203 j = self.rng.randint(0, high=n)
204 ind1 = np.where(tags == i)
205 ind2 = np.where(tags == j)
206 cop1 = np.mean(pos[ind1], axis=0)
207 cop2 = np.mean(pos[ind2], axis=0)
208 pos[ind1] += cop2 - cop1
209 pos[ind2] += cop1 - cop2
211 top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
212 if self.blmin is None:
213 too_close = False
214 else:
215 too_close = atoms_too_close(
216 top, self.blmin, use_tags=self.use_tags)
217 if not too_close and self.test_dist_to_slab:
218 too_close = atoms_too_close_two_sets(top, slab, self.blmin)
220 if count == maxcount:
221 return None
223 mutant = slab + top
224 return mutant
227class MirrorMutation(OffspringCreator):
228 """A mirror mutation, as described in
229 TO BE PUBLISHED.
231 This mutation mirrors half of the cluster in a
232 randomly oriented cutting plane discarding the other half.
234 Parameters:
236 blmin: Dictionary defining the minimum allowed
237 distance between atoms.
239 n_top: Number of atoms the GA optimizes.
241 reflect: Defines if the mirrored half is also reflected
242 perpendicular to the mirroring plane.
244 rng: Random number generator
245 By default numpy.random.
246 """
248 def __init__(self, blmin, n_top, reflect=False, rng=np.random,
249 verbose=False):
250 OffspringCreator.__init__(self, verbose, rng=rng)
251 self.blmin = blmin
252 self.n_top = n_top
253 self.reflect = reflect
255 self.descriptor = 'MirrorMutation'
256 self.min_inputs = 1
258 def get_new_individual(self, parents):
259 f = parents[0]
261 indi = self.mutate(f)
262 if indi is None:
263 return indi, 'mutation: mirror'
265 indi = self.initialize_individual(f, indi)
266 indi.info['data']['parents'] = [f.info['confid']]
268 return self.finalize_individual(indi), 'mutation: mirror'
270 def mutate(self, atoms):
271 """ Do the mutation of the atoms input. """
272 reflect = self.reflect
273 tc = True
274 slab = atoms[0:len(atoms) - self.n_top]
275 top = atoms[len(atoms) - self.n_top: len(atoms)]
276 num = top.numbers
277 unique_types = list(set(num))
278 nu = {u: sum(num == u) for u in unique_types}
279 n_tries = 1000
280 counter = 0
281 changed = False
283 while tc and counter < n_tries:
284 counter += 1
285 cand = top.copy()
286 pos = cand.get_positions()
288 cm = np.average(top.get_positions(), axis=0)
290 # first select a randomly oriented cutting plane
291 theta = pi * self.rng.random()
292 phi = 2. * pi * self.rng.random()
293 n = (cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta))
294 n = np.array(n)
296 # Calculate all atoms signed distance to the cutting plane
297 D = []
298 for (i, p) in enumerate(pos):
299 d = np.dot(p - cm, n)
300 D.append((i, d))
302 # Sort the atoms by their signed distance
303 D.sort(key=lambda x: x[1])
304 nu_taken = {}
306 # Select half of the atoms needed for a full cluster
307 p_use = []
308 n_use = []
309 for (i, d) in D:
310 if num[i] not in nu_taken.keys():
311 nu_taken[num[i]] = 0
312 if nu_taken[num[i]] < nu[num[i]] / 2.:
313 p_use.append(pos[i])
314 n_use.append(num[i])
315 nu_taken[num[i]] += 1
317 # calculate the mirrored position and add these.
318 pn = []
319 for p in p_use:
320 pt = p - 2. * np.dot(p - cm, n) * n
321 if reflect:
322 pt = -pt + 2 * cm + 2 * n * np.dot(pt - cm, n)
323 pn.append(pt)
325 n_use.extend(n_use)
326 p_use.extend(pn)
328 # In the case of an uneven number of
329 # atoms we need to add one extra
330 for n in nu:
331 if nu[n] % 2 == 0:
332 continue
333 while sum(n_use == n) > nu[n]:
334 for i in range(int(len(n_use) / 2), len(n_use)):
335 if n_use[i] == n:
336 del p_use[i]
337 del n_use[i]
338 break
339 assert sum(n_use == n) == nu[n]
341 # Make sure we have the correct number of atoms
342 # and rearrange the atoms so they are in the right order
343 for i in range(len(n_use)):
344 if num[i] == n_use[i]:
345 continue
346 for j in range(i + 1, len(n_use)):
347 if n_use[j] == num[i]:
348 tn = n_use[i]
349 tp = p_use[i]
350 n_use[i] = n_use[j]
351 p_use[i] = p_use[j]
352 p_use[j] = tp
353 n_use[j] = tn
355 # Finally we check that nothing is too close in the end product.
356 cand = Atoms(num, p_use, cell=slab.get_cell(), pbc=slab.get_pbc())
358 tc = atoms_too_close(cand, self.blmin)
359 if tc:
360 continue
361 tc = atoms_too_close_two_sets(slab, cand, self.blmin)
363 if not changed and counter > n_tries // 2:
364 reflect = not reflect
365 changed = True
367 tot = slab + cand
369 if counter == n_tries:
370 return None
371 return tot
374class StrainMutation(OffspringCreator):
375 """ Mutates a candidate by applying a randomly generated strain.
377 For more information, see also:
379 * :doi:`Glass, Oganov, Hansen, Comp. Phys. Comm. 175 (2006) 713-720
380 <10.1016/j.cpc.2006.07.020>`
382 * :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
383 <10.1016/j.cpc.2010.07.048>`
385 After initialization of the mutation, a scaling volume
386 (to which each mutated structure is scaled before checking the
387 constraints) is typically generated from the population,
388 which is then also occasionally updated in the course of the
389 GA run.
391 Parameters:
393 blmin: dict
394 The closest allowed interatomic distances on the form:
395 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
397 cellbounds: ase.ga.utilities.CellBounds instance
398 Describes limits on the cell shape, see
399 :class:`~ase.ga.utilities.CellBounds`.
401 stddev: float
402 Standard deviation used in the generation of the
403 strain matrix elements.
405 number_of_variable_cell_vectors: int (default 3)
406 The number of variable cell vectors (1, 2 or 3).
407 To keep things simple, it is the 'first' vectors which
408 will be treated as variable, i.e. the 'a' vector in the
409 univariate case, the 'a' and 'b' vectors in the bivariate
410 case, etc.
412 use_tags: boolean
413 Whether to use the atomic tags to preserve molecular identity.
415 rng: Random number generator
416 By default numpy.random.
417 """
419 def __init__(self, blmin, cellbounds=None, stddev=0.7,
420 number_of_variable_cell_vectors=3, use_tags=False,
421 rng=np.random, verbose=False):
422 OffspringCreator.__init__(self, verbose, rng=rng)
423 self.blmin = blmin
424 self.cellbounds = cellbounds
425 self.stddev = stddev
426 self.number_of_variable_cell_vectors = number_of_variable_cell_vectors
427 self.use_tags = use_tags
429 self.scaling_volume = None
430 self.descriptor = 'StrainMutation'
431 self.min_inputs = 1
433 def update_scaling_volume(self, population, w_adapt=0.5, n_adapt=0):
434 """Function to initialize or update the scaling volume in a GA run.
436 w_adapt: weight of the new vs the old scaling volume
438 n_adapt: number of best candidates in the population that
439 are used to calculate the new scaling volume
440 """
441 if not n_adapt:
442 # if not set, take best 20% of the population
443 n_adapt = int(np.ceil(0.2 * len(population)))
444 v_new = np.mean([a.get_volume() for a in population[:n_adapt]])
446 if not self.scaling_volume:
447 self.scaling_volume = v_new
448 else:
449 volumes = [self.scaling_volume, v_new]
450 weights = [1 - w_adapt, w_adapt]
451 self.scaling_volume = np.average(volumes, weights=weights)
453 def get_new_individual(self, parents):
454 f = parents[0]
456 indi = self.mutate(f)
457 if indi is None:
458 return indi, 'mutation: strain'
460 indi = self.initialize_individual(f, indi)
461 indi.info['data']['parents'] = [f.info['confid']]
463 return self.finalize_individual(indi), 'mutation: strain'
465 def mutate(self, atoms):
466 """ Does the actual mutation. """
467 cell_ref = atoms.get_cell()
468 pos_ref = atoms.get_positions()
470 if self.scaling_volume is None:
471 # The scaling_volume has not been set (yet),
472 # so we give it the same volume as the parent
473 vol_ref = atoms.get_volume()
474 else:
475 vol_ref = self.scaling_volume
477 if self.use_tags:
478 tags = atoms.get_tags()
479 gather_atoms_by_tag(atoms)
480 pos = atoms.get_positions()
482 mutant = atoms.copy()
484 count = 0
485 too_close = True
486 maxcount = 1000
487 while too_close and count < maxcount:
488 count += 1
490 # generating the strain matrix:
491 strain = np.identity(3)
492 for i in range(self.number_of_variable_cell_vectors):
493 for j in range(i + 1):
494 r = self.rng.normal(loc=0., scale=self.stddev)
495 if i == j:
496 strain[i, j] += r
497 else:
498 epsilon = 0.5 * r
499 strain[i, j] += epsilon
500 strain[j, i] += epsilon
502 # applying the strain:
503 cell_new = np.dot(strain, cell_ref)
505 # convert the submatrix with the variable cell vectors
506 # to a lower triangular form
507 cell_new = calc_rotated_cell(cell_new)
508 for i in range(self.number_of_variable_cell_vectors, 3):
509 cell_new[i] = cell_ref[i]
511 cell_new = Cell(cell_new)
513 # volume scaling:
514 if self.number_of_variable_cell_vectors > 0:
515 scaling = vol_ref / cell_new.volume
516 scaling **= 1. / self.number_of_variable_cell_vectors
517 cell_new[:self.number_of_variable_cell_vectors] *= scaling
519 # check cell dimensions:
520 if not self.cellbounds.is_within_bounds(cell_new):
521 continue
523 # ensure non-variable cell vectors are indeed unchanged
524 for i in range(self.number_of_variable_cell_vectors, 3):
525 assert np.allclose(cell_new[i], cell_ref[i])
527 # check that the volume is correct
528 assert np.isclose(vol_ref, cell_new.volume)
530 # apply the new unit cell and scale
531 # the atomic positions accordingly
532 mutant.set_cell(cell_ref, scale_atoms=False)
534 if self.use_tags:
535 transfo = np.linalg.solve(cell_ref, cell_new)
536 for tag in np.unique(tags):
537 select = np.where(tags == tag)
538 cop = np.mean(pos[select], axis=0)
539 disp = np.dot(cop, transfo) - cop
540 mutant.positions[select] += disp
541 else:
542 mutant.set_positions(pos_ref)
544 mutant.set_cell(cell_new, scale_atoms=not self.use_tags)
545 mutant.wrap()
547 # check the interatomic distances
548 too_close = atoms_too_close(mutant, self.blmin,
549 use_tags=self.use_tags)
551 if count == maxcount:
552 mutant = None
554 return mutant
557class PermuStrainMutation(CombinationMutation):
558 """Combination of PermutationMutation and StrainMutation.
560 For more information, see also:
562 * :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
563 <10.1016/j.cpc.2010.07.048>`
565 Parameters:
567 permutationmutation: OffspringCreator instance
568 A mutation that permutes atom types.
570 strainmutation: OffspringCreator instance
571 A mutation that mutates by straining.
572 """
574 def __init__(self, permutationmutation, strainmutation, verbose=False):
575 super().__init__(permutationmutation,
576 strainmutation,
577 verbose=verbose)
578 self.descriptor = 'permustrain'
581class RotationalMutation(OffspringCreator):
582 """Mutates a candidate by applying random rotations
583 to multi-atom moieties in the structure (atoms with
584 the same tag are considered part of one such moiety).
586 Only performs whole-molecule rotations, no internal
587 rotations.
589 For more information, see also:
591 * `Zhu Q., Oganov A.R., Glass C.W., Stokes H.T,
592 Acta Cryst. (2012), B68, 215-226.`__
594 __ https://dx.doi.org/10.1107/S0108768112017466
596 Parameters:
598 blmin: dict
599 The closest allowed interatomic distances on the form:
600 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
602 n_top: int or None
603 The number of atoms to optimize (None = include all).
605 fraction: float
606 Fraction of the moieties to be rotated.
608 tags: None or list of integers
609 Specifies, respectively, whether all moieties or only those
610 with matching tags are eligible for rotation.
612 min_angle: float
613 Minimal angle (in radians) for each rotation;
614 should lie in the interval [0, pi].
616 test_dist_to_slab: boolean
617 Whether also the distances to the slab
618 should be checked to satisfy the blmin.
620 rng: Random number generator
621 By default numpy.random.
622 """
624 def __init__(self, blmin, n_top=None, fraction=0.33, tags=None,
625 min_angle=1.57, test_dist_to_slab=True, rng=np.random,
626 verbose=False):
627 OffspringCreator.__init__(self, verbose, rng=rng)
628 self.blmin = blmin
629 self.n_top = n_top
630 self.fraction = fraction
631 self.tags = tags
632 self.min_angle = min_angle
633 self.test_dist_to_slab = test_dist_to_slab
634 self.descriptor = 'RotationalMutation'
635 self.min_inputs = 1
637 def get_new_individual(self, parents):
638 f = parents[0]
640 indi = self.mutate(f)
641 if indi is None:
642 return indi, 'mutation: rotational'
644 indi = self.initialize_individual(f, indi)
645 indi.info['data']['parents'] = [f.info['confid']]
647 return self.finalize_individual(indi), 'mutation: rotational'
649 def mutate(self, atoms):
650 """Does the actual mutation."""
651 N = len(atoms) if self.n_top is None else self.n_top
652 slab = atoms[:len(atoms) - N]
653 atoms = atoms[-N:]
655 mutant = atoms.copy()
656 gather_atoms_by_tag(mutant)
657 pos = mutant.get_positions()
658 tags = mutant.get_tags()
659 eligible_tags = tags if self.tags is None else self.tags
661 indices = {}
662 for tag in np.unique(tags):
663 hits = np.where(tags == tag)[0]
664 if len(hits) > 1 and tag in eligible_tags:
665 indices[tag] = hits
667 n_rot = int(np.ceil(len(indices) * self.fraction))
668 chosen_tags = self.rng.choice(list(indices.keys()), size=n_rot,
669 replace=False)
671 too_close = True
672 count = 0
673 maxcount = 10000
674 while too_close and count < maxcount:
675 newpos = np.copy(pos)
676 for tag in chosen_tags:
677 p = np.copy(newpos[indices[tag]])
678 cop = np.mean(p, axis=0)
680 if len(p) == 2:
681 line = (p[1] - p[0]) / np.linalg.norm(p[1] - p[0])
682 while True:
683 axis = self.rng.random(3)
684 axis /= np.linalg.norm(axis)
685 a = np.arccos(np.dot(axis, line))
686 if np.pi / 4 < a < np.pi * 3 / 4:
687 break
688 else:
689 axis = self.rng.random(3)
690 axis /= np.linalg.norm(axis)
692 angle = self.min_angle
693 angle += 2 * (np.pi - self.min_angle) * self.rng.random()
695 m = get_rotation_matrix(axis, angle)
696 newpos[indices[tag]] = np.dot(m, (p - cop).T).T + cop
698 mutant.set_positions(newpos)
699 mutant.wrap()
700 too_close = atoms_too_close(mutant, self.blmin, use_tags=True)
701 count += 1
703 if not too_close and self.test_dist_to_slab:
704 too_close = atoms_too_close_two_sets(slab, mutant, self.blmin)
706 if count == maxcount:
707 mutant = None
708 else:
709 mutant = slab + mutant
711 return mutant
714class RattleRotationalMutation(CombinationMutation):
715 """Combination of RattleMutation and RotationalMutation.
717 Parameters:
719 rattlemutation: OffspringCreator instance
720 A mutation that rattles atoms.
722 rotationalmutation: OffspringCreator instance
723 A mutation that rotates moieties.
724 """
726 def __init__(self, rattlemutation, rotationalmutation, verbose=False):
727 super().__init__(rattlemutation,
728 rotationalmutation,
729 verbose=verbose)
730 self.descriptor = 'rattlerotational'