Coverage for /builds/kinetik161/ase/ase/ga/standardmutations.py: 73.67%
376 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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 (atoms_too_close, atoms_too_close_two_sets,
11 gather_atoms_by_tag, get_rotation_matrix)
14class RattleMutation(OffspringCreator):
15 """An implementation of the rattle mutation as described in:
17 R.L. Johnston Dalton Transactions, Vol. 22,
18 No. 22. (2003), pp. 4193-4207
20 Parameters:
22 blmin: Dictionary defining the minimum distance between atoms
23 after the rattle.
25 n_top: Number of atoms optimized by the GA.
27 rattle_strength: Strength with which the atoms are moved.
29 rattle_prop: The probability with which each atom is rattled.
31 test_dist_to_slab: whether to also make sure that the distances
32 between the atoms and the slab satisfy the blmin.
34 use_tags: if True, the atomic tags will be used to preserve
35 molecular identity. Same-tag atoms will then be
36 displaced collectively, so that the internal
37 geometry is preserved.
39 rng: Random number generator
40 By default numpy.random.
41 """
43 def __init__(self, blmin, n_top, rattle_strength=0.8,
44 rattle_prop=0.4, test_dist_to_slab=True, use_tags=False,
45 verbose=False, rng=np.random):
46 OffspringCreator.__init__(self, verbose, rng=rng)
47 self.blmin = blmin
48 self.n_top = n_top
49 self.rattle_strength = rattle_strength
50 self.rattle_prop = rattle_prop
51 self.test_dist_to_slab = test_dist_to_slab
52 self.use_tags = use_tags
54 self.descriptor = 'RattleMutation'
55 self.min_inputs = 1
57 def get_new_individual(self, parents):
58 f = parents[0]
60 indi = self.mutate(f)
61 if indi is None:
62 return indi, 'mutation: rattle'
64 indi = self.initialize_individual(f, indi)
65 indi.info['data']['parents'] = [f.info['confid']]
67 return self.finalize_individual(indi), 'mutation: rattle'
69 def mutate(self, atoms):
70 """Does the actual mutation."""
71 N = len(atoms) if self.n_top is None else self.n_top
72 slab = atoms[:len(atoms) - N]
73 atoms = atoms[-N:]
74 tags = atoms.get_tags() if self.use_tags else np.arange(N)
75 pos_ref = atoms.get_positions()
76 num = atoms.get_atomic_numbers()
77 cell = atoms.get_cell()
78 pbc = atoms.get_pbc()
79 st = 2. * self.rattle_strength
81 count = 0
82 maxcount = 1000
83 too_close = True
84 while too_close and count < maxcount:
85 count += 1
86 pos = pos_ref.copy()
87 ok = False
88 for tag in np.unique(tags):
89 select = np.where(tags == tag)
90 if self.rng.random() < self.rattle_prop:
91 ok = True
92 r = self.rng.random(3)
93 pos[select] += st * (r - 0.5)
95 if not ok:
96 # Nothing got rattled
97 continue
99 top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
100 too_close = atoms_too_close(
101 top, self.blmin, use_tags=self.use_tags)
102 if not too_close and self.test_dist_to_slab:
103 too_close = atoms_too_close_two_sets(top, slab, self.blmin)
105 if count == maxcount:
106 return None
108 mutant = slab + top
109 return mutant
112class PermutationMutation(OffspringCreator):
113 """Mutation that permutes a percentage of the atom types in the cluster.
115 Parameters:
117 n_top: Number of atoms optimized by the GA.
119 probability: The probability with which an atom is permuted.
121 test_dist_to_slab: whether to also make sure that the distances
122 between the atoms and the slab satisfy the blmin.
124 use_tags: if True, the atomic tags will be used to preserve
125 molecular identity. Permutations will then happen
126 at the molecular level, i.e. swapping the center-of-
127 positions of two moieties while preserving their
128 internal geometries.
130 blmin: Dictionary defining the minimum distance between atoms
131 after the permutation. If equal to None (the default),
132 no such check is performed.
134 rng: Random number generator
135 By default numpy.random.
136 """
138 def __init__(self, n_top, probability=0.33, test_dist_to_slab=True,
139 use_tags=False, blmin=None, rng=np.random, verbose=False):
140 OffspringCreator.__init__(self, verbose, rng=rng)
141 self.n_top = n_top
142 self.probability = probability
143 self.test_dist_to_slab = test_dist_to_slab
144 self.use_tags = use_tags
145 self.blmin = blmin
147 self.descriptor = 'PermutationMutation'
148 self.min_inputs = 1
150 def get_new_individual(self, parents):
151 f = parents[0]
153 indi = self.mutate(f)
154 if indi is None:
155 return indi, 'mutation: permutation'
157 indi = self.initialize_individual(f, indi)
158 indi.info['data']['parents'] = [f.info['confid']]
160 return self.finalize_individual(indi), 'mutation: permutation'
162 def mutate(self, atoms):
163 """Does the actual mutation."""
164 N = len(atoms) if self.n_top is None else self.n_top
165 slab = atoms[:len(atoms) - N]
166 atoms = atoms[-N:]
167 if self.use_tags:
168 gather_atoms_by_tag(atoms)
169 tags = atoms.get_tags() if self.use_tags else np.arange(N)
170 pos_ref = atoms.get_positions()
171 num = atoms.get_atomic_numbers()
172 cell = atoms.get_cell()
173 pbc = atoms.get_pbc()
174 symbols = atoms.get_chemical_symbols()
176 unique_tags = np.unique(tags)
177 n = len(unique_tags)
178 swaps = int(np.ceil(n * self.probability / 2.))
180 sym = []
181 for tag in unique_tags:
182 indices = np.where(tags == tag)[0]
183 s = ''.join([symbols[j] for j in indices])
184 sym.append(s)
186 assert len(np.unique(sym)) > 1, \
187 'Permutations with one atom (or molecule) type is not valid'
189 count = 0
190 maxcount = 1000
191 too_close = True
192 while too_close and count < maxcount:
193 count += 1
194 pos = pos_ref.copy()
195 for _ in range(swaps):
196 i = j = 0
197 while sym[i] == sym[j]:
198 i = self.rng.randint(0, high=n)
199 j = self.rng.randint(0, high=n)
200 ind1 = np.where(tags == i)
201 ind2 = np.where(tags == j)
202 cop1 = np.mean(pos[ind1], axis=0)
203 cop2 = np.mean(pos[ind2], axis=0)
204 pos[ind1] += cop2 - cop1
205 pos[ind2] += cop1 - cop2
207 top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
208 if self.blmin is None:
209 too_close = False
210 else:
211 too_close = atoms_too_close(
212 top, self.blmin, use_tags=self.use_tags)
213 if not too_close and self.test_dist_to_slab:
214 too_close = atoms_too_close_two_sets(top, slab, self.blmin)
216 if count == maxcount:
217 return None
219 mutant = slab + top
220 return mutant
223class MirrorMutation(OffspringCreator):
224 """A mirror mutation, as described in
225 TO BE PUBLISHED.
227 This mutation mirrors half of the cluster in a
228 randomly oriented cutting plane discarding the other half.
230 Parameters:
232 blmin: Dictionary defining the minimum allowed
233 distance between atoms.
235 n_top: Number of atoms the GA optimizes.
237 reflect: Defines if the mirrored half is also reflected
238 perpendicular to the mirroring plane.
240 rng: Random number generator
241 By default numpy.random.
242 """
244 def __init__(self, blmin, n_top, reflect=False, rng=np.random,
245 verbose=False):
246 OffspringCreator.__init__(self, verbose, rng=rng)
247 self.blmin = blmin
248 self.n_top = n_top
249 self.reflect = reflect
251 self.descriptor = 'MirrorMutation'
252 self.min_inputs = 1
254 def get_new_individual(self, parents):
255 f = parents[0]
257 indi = self.mutate(f)
258 if indi is None:
259 return indi, 'mutation: mirror'
261 indi = self.initialize_individual(f, indi)
262 indi.info['data']['parents'] = [f.info['confid']]
264 return self.finalize_individual(indi), 'mutation: mirror'
266 def mutate(self, atoms):
267 """ Do the mutation of the atoms input. """
268 reflect = self.reflect
269 tc = True
270 slab = atoms[0:len(atoms) - self.n_top]
271 top = atoms[len(atoms) - self.n_top: len(atoms)]
272 num = top.numbers
273 unique_types = list(set(num))
274 nu = {}
275 for u in unique_types:
276 nu[u] = sum(num == u)
278 n_tries = 1000
279 counter = 0
280 changed = False
282 while tc and counter < n_tries:
283 counter += 1
284 cand = top.copy()
285 pos = cand.get_positions()
287 cm = np.average(top.get_positions(), axis=0)
289 # first select a randomly oriented cutting plane
290 theta = pi * self.rng.random()
291 phi = 2. * pi * self.rng.random()
292 n = (cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta))
293 n = np.array(n)
295 # Calculate all atoms signed distance to the cutting plane
296 D = []
297 for (i, p) in enumerate(pos):
298 d = np.dot(p - cm, n)
299 D.append((i, d))
301 # Sort the atoms by their signed distance
302 D.sort(key=lambda x: x[1])
303 nu_taken = {}
305 # Select half of the atoms needed for a full cluster
306 p_use = []
307 n_use = []
308 for (i, d) in D:
309 if num[i] not in nu_taken.keys():
310 nu_taken[num[i]] = 0
311 if nu_taken[num[i]] < nu[num[i]] / 2.:
312 p_use.append(pos[i])
313 n_use.append(num[i])
314 nu_taken[num[i]] += 1
316 # calculate the mirrored position and add these.
317 pn = []
318 for p in p_use:
319 pt = p - 2. * np.dot(p - cm, n) * n
320 if reflect:
321 pt = -pt + 2 * cm + 2 * n * np.dot(pt - cm, n)
322 pn.append(pt)
324 n_use.extend(n_use)
325 p_use.extend(pn)
327 # In the case of an uneven number of
328 # atoms we need to add one extra
329 for n in nu:
330 if nu[n] % 2 == 0:
331 continue
332 while sum(n_use == n) > nu[n]:
333 for i in range(int(len(n_use) / 2), len(n_use)):
334 if n_use[i] == n:
335 del p_use[i]
336 del n_use[i]
337 break
338 assert sum(n_use == n) == nu[n]
340 # Make sure we have the correct number of atoms
341 # and rearrange the atoms so they are in the right order
342 for i in range(len(n_use)):
343 if num[i] == n_use[i]:
344 continue
345 for j in range(i + 1, len(n_use)):
346 if n_use[j] == num[i]:
347 tn = n_use[i]
348 tp = p_use[i]
349 n_use[i] = n_use[j]
350 p_use[i] = p_use[j]
351 p_use[j] = tp
352 n_use[j] = tn
354 # Finally we check that nothing is too close in the end product.
355 cand = Atoms(num, p_use, cell=slab.get_cell(), pbc=slab.get_pbc())
357 tc = atoms_too_close(cand, self.blmin)
358 if tc:
359 continue
360 tc = atoms_too_close_two_sets(slab, cand, self.blmin)
362 if not changed and counter > n_tries // 2:
363 reflect = not reflect
364 changed = True
366 tot = slab + cand
368 if counter == n_tries:
369 return None
370 return tot
373class StrainMutation(OffspringCreator):
374 """ Mutates a candidate by applying a randomly generated strain.
376 For more information, see also:
378 * :doi:`Glass, Oganov, Hansen, Comp. Phys. Comm. 175 (2006) 713-720
379 <10.1016/j.cpc.2006.07.020>`
381 * :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
382 <10.1016/j.cpc.2010.07.048>`
384 After initialization of the mutation, a scaling volume
385 (to which each mutated structure is scaled before checking the
386 constraints) is typically generated from the population,
387 which is then also occasionally updated in the course of the
388 GA run.
390 Parameters:
392 blmin: dict
393 The closest allowed interatomic distances on the form:
394 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
396 cellbounds: ase.ga.utilities.CellBounds instance
397 Describes limits on the cell shape, see
398 :class:`~ase.ga.utilities.CellBounds`.
400 stddev: float
401 Standard deviation used in the generation of the
402 strain matrix elements.
404 number_of_variable_cell_vectors: int (default 3)
405 The number of variable cell vectors (1, 2 or 3).
406 To keep things simple, it is the 'first' vectors which
407 will be treated as variable, i.e. the 'a' vector in the
408 univariate case, the 'a' and 'b' vectors in the bivariate
409 case, etc.
411 use_tags: boolean
412 Whether to use the atomic tags to preserve molecular identity.
414 rng: Random number generator
415 By default numpy.random.
416 """
418 def __init__(self, blmin, cellbounds=None, stddev=0.7,
419 number_of_variable_cell_vectors=3, use_tags=False,
420 rng=np.random, verbose=False):
421 OffspringCreator.__init__(self, verbose, rng=rng)
422 self.blmin = blmin
423 self.cellbounds = cellbounds
424 self.stddev = stddev
425 self.number_of_variable_cell_vectors = number_of_variable_cell_vectors
426 self.use_tags = use_tags
428 self.scaling_volume = None
429 self.descriptor = 'StrainMutation'
430 self.min_inputs = 1
432 def update_scaling_volume(self, population, w_adapt=0.5, n_adapt=0):
433 """Function to initialize or update the scaling volume in a GA run.
435 w_adapt: weight of the new vs the old scaling volume
437 n_adapt: number of best candidates in the population that
438 are used to calculate the new scaling volume
439 """
440 if not n_adapt:
441 # if not set, take best 20% of the population
442 n_adapt = int(np.ceil(0.2 * len(population)))
443 v_new = np.mean([a.get_volume() for a in population[:n_adapt]])
445 if not self.scaling_volume:
446 self.scaling_volume = v_new
447 else:
448 volumes = [self.scaling_volume, v_new]
449 weights = [1 - w_adapt, w_adapt]
450 self.scaling_volume = np.average(volumes, weights=weights)
452 def get_new_individual(self, parents):
453 f = parents[0]
455 indi = self.mutate(f)
456 if indi is None:
457 return indi, 'mutation: strain'
459 indi = self.initialize_individual(f, indi)
460 indi.info['data']['parents'] = [f.info['confid']]
462 return self.finalize_individual(indi), 'mutation: strain'
464 def mutate(self, atoms):
465 """ Does the actual mutation. """
466 cell_ref = atoms.get_cell()
467 pos_ref = atoms.get_positions()
469 if self.scaling_volume is None:
470 # The scaling_volume has not been set (yet),
471 # so we give it the same volume as the parent
472 vol_ref = atoms.get_volume()
473 else:
474 vol_ref = self.scaling_volume
476 if self.use_tags:
477 tags = atoms.get_tags()
478 gather_atoms_by_tag(atoms)
479 pos = atoms.get_positions()
481 mutant = atoms.copy()
483 count = 0
484 too_close = True
485 maxcount = 1000
486 while too_close and count < maxcount:
487 count += 1
489 # generating the strain matrix:
490 strain = np.identity(3)
491 for i in range(self.number_of_variable_cell_vectors):
492 for j in range(i + 1):
493 r = self.rng.normal(loc=0., scale=self.stddev)
494 if i == j:
495 strain[i, j] += r
496 else:
497 epsilon = 0.5 * r
498 strain[i, j] += epsilon
499 strain[j, i] += epsilon
501 # applying the strain:
502 cell_new = np.dot(strain, cell_ref)
504 # convert the submatrix with the variable cell vectors
505 # to a lower triangular form
506 cell_new = calc_rotated_cell(cell_new)
507 for i in range(self.number_of_variable_cell_vectors, 3):
508 cell_new[i] = cell_ref[i]
510 cell_new = Cell(cell_new)
512 # volume scaling:
513 if self.number_of_variable_cell_vectors > 0:
514 scaling = vol_ref / cell_new.volume
515 scaling **= 1. / self.number_of_variable_cell_vectors
516 cell_new[:self.number_of_variable_cell_vectors] *= scaling
518 # check cell dimensions:
519 if not self.cellbounds.is_within_bounds(cell_new):
520 continue
522 # ensure non-variable cell vectors are indeed unchanged
523 for i in range(self.number_of_variable_cell_vectors, 3):
524 assert np.allclose(cell_new[i], cell_ref[i])
526 # check that the volume is correct
527 assert np.isclose(vol_ref, cell_new.volume)
529 # apply the new unit cell and scale
530 # the atomic positions accordingly
531 mutant.set_cell(cell_ref, scale_atoms=False)
533 if self.use_tags:
534 transfo = np.linalg.solve(cell_ref, cell_new)
535 for tag in np.unique(tags):
536 select = np.where(tags == tag)
537 cop = np.mean(pos[select], axis=0)
538 disp = np.dot(cop, transfo) - cop
539 mutant.positions[select] += disp
540 else:
541 mutant.set_positions(pos_ref)
543 mutant.set_cell(cell_new, scale_atoms=not self.use_tags)
544 mutant.wrap()
546 # check the interatomic distances
547 too_close = atoms_too_close(mutant, self.blmin,
548 use_tags=self.use_tags)
550 if count == maxcount:
551 mutant = None
553 return mutant
556class PermuStrainMutation(CombinationMutation):
557 """Combination of PermutationMutation and StrainMutation.
559 For more information, see also:
561 * :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
562 <10.1016/j.cpc.2010.07.048>`
564 Parameters:
566 permutationmutation: OffspringCreator instance
567 A mutation that permutes atom types.
569 strainmutation: OffspringCreator instance
570 A mutation that mutates by straining.
571 """
573 def __init__(self, permutationmutation, strainmutation, verbose=False):
574 super().__init__(permutationmutation,
575 strainmutation,
576 verbose=verbose)
577 self.descriptor = 'permustrain'
580class RotationalMutation(OffspringCreator):
581 """Mutates a candidate by applying random rotations
582 to multi-atom moieties in the structure (atoms with
583 the same tag are considered part of one such moiety).
585 Only performs whole-molecule rotations, no internal
586 rotations.
588 For more information, see also:
590 * `Zhu Q., Oganov A.R., Glass C.W., Stokes H.T,
591 Acta Cryst. (2012), B68, 215-226.`__
593 __ https://dx.doi.org/10.1107/S0108768112017466
595 Parameters:
597 blmin: dict
598 The closest allowed interatomic distances on the form:
599 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
601 n_top: int or None
602 The number of atoms to optimize (None = include all).
604 fraction: float
605 Fraction of the moieties to be rotated.
607 tags: None or list of integers
608 Specifies, respectively, whether all moieties or only those
609 with matching tags are eligible for rotation.
611 min_angle: float
612 Minimal angle (in radians) for each rotation;
613 should lie in the interval [0, pi].
615 test_dist_to_slab: boolean
616 Whether also the distances to the slab
617 should be checked to satisfy the blmin.
619 rng: Random number generator
620 By default numpy.random.
621 """
623 def __init__(self, blmin, n_top=None, fraction=0.33, tags=None,
624 min_angle=1.57, test_dist_to_slab=True, rng=np.random,
625 verbose=False):
626 OffspringCreator.__init__(self, verbose, rng=rng)
627 self.blmin = blmin
628 self.n_top = n_top
629 self.fraction = fraction
630 self.tags = tags
631 self.min_angle = min_angle
632 self.test_dist_to_slab = test_dist_to_slab
633 self.descriptor = 'RotationalMutation'
634 self.min_inputs = 1
636 def get_new_individual(self, parents):
637 f = parents[0]
639 indi = self.mutate(f)
640 if indi is None:
641 return indi, 'mutation: rotational'
643 indi = self.initialize_individual(f, indi)
644 indi.info['data']['parents'] = [f.info['confid']]
646 return self.finalize_individual(indi), 'mutation: rotational'
648 def mutate(self, atoms):
649 """Does the actual mutation."""
650 N = len(atoms) if self.n_top is None else self.n_top
651 slab = atoms[:len(atoms) - N]
652 atoms = atoms[-N:]
654 mutant = atoms.copy()
655 gather_atoms_by_tag(mutant)
656 pos = mutant.get_positions()
657 tags = mutant.get_tags()
658 eligible_tags = tags if self.tags is None else self.tags
660 indices = {}
661 for tag in np.unique(tags):
662 hits = np.where(tags == tag)[0]
663 if len(hits) > 1 and tag in eligible_tags:
664 indices[tag] = hits
666 n_rot = int(np.ceil(len(indices) * self.fraction))
667 chosen_tags = self.rng.choice(list(indices.keys()), size=n_rot,
668 replace=False)
670 too_close = True
671 count = 0
672 maxcount = 10000
673 while too_close and count < maxcount:
674 newpos = np.copy(pos)
675 for tag in chosen_tags:
676 p = np.copy(newpos[indices[tag]])
677 cop = np.mean(p, axis=0)
679 if len(p) == 2:
680 line = (p[1] - p[0]) / np.linalg.norm(p[1] - p[0])
681 while True:
682 axis = self.rng.random(3)
683 axis /= np.linalg.norm(axis)
684 a = np.arccos(np.dot(axis, line))
685 if np.pi / 4 < a < np.pi * 3 / 4:
686 break
687 else:
688 axis = self.rng.random(3)
689 axis /= np.linalg.norm(axis)
691 angle = self.min_angle
692 angle += 2 * (np.pi - self.min_angle) * self.rng.random()
694 m = get_rotation_matrix(axis, angle)
695 newpos[indices[tag]] = np.dot(m, (p - cop).T).T + cop
697 mutant.set_positions(newpos)
698 mutant.wrap()
699 too_close = atoms_too_close(mutant, self.blmin, use_tags=True)
700 count += 1
702 if not too_close and self.test_dist_to_slab:
703 too_close = atoms_too_close_two_sets(slab, mutant, self.blmin)
705 if count == maxcount:
706 mutant = None
707 else:
708 mutant = slab + mutant
710 return mutant
713class RattleRotationalMutation(CombinationMutation):
714 """Combination of RattleMutation and RotationalMutation.
716 Parameters:
718 rattlemutation: OffspringCreator instance
719 A mutation that rattles atoms.
721 rotationalmutation: OffspringCreator instance
722 A mutation that rotates moieties.
723 """
725 def __init__(self, rattlemutation, rotationalmutation, verbose=False):
726 super().__init__(rattlemutation,
727 rotationalmutation,
728 verbose=verbose)
729 self.descriptor = 'rattlerotational'