Coverage for /builds/hweiske/ase/ase/ga/standardmutations.py: 74.06%
374 statements
« prev ^ index » next coverage.py v7.2.7, created at 2024-04-22 11:22 +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 = {u: sum(num == u) for u in unique_types}
275 n_tries = 1000
276 counter = 0
277 changed = False
279 while tc and counter < n_tries:
280 counter += 1
281 cand = top.copy()
282 pos = cand.get_positions()
284 cm = np.average(top.get_positions(), axis=0)
286 # first select a randomly oriented cutting plane
287 theta = pi * self.rng.random()
288 phi = 2. * pi * self.rng.random()
289 n = (cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta))
290 n = np.array(n)
292 # Calculate all atoms signed distance to the cutting plane
293 D = []
294 for (i, p) in enumerate(pos):
295 d = np.dot(p - cm, n)
296 D.append((i, d))
298 # Sort the atoms by their signed distance
299 D.sort(key=lambda x: x[1])
300 nu_taken = {}
302 # Select half of the atoms needed for a full cluster
303 p_use = []
304 n_use = []
305 for (i, d) in D:
306 if num[i] not in nu_taken.keys():
307 nu_taken[num[i]] = 0
308 if nu_taken[num[i]] < nu[num[i]] / 2.:
309 p_use.append(pos[i])
310 n_use.append(num[i])
311 nu_taken[num[i]] += 1
313 # calculate the mirrored position and add these.
314 pn = []
315 for p in p_use:
316 pt = p - 2. * np.dot(p - cm, n) * n
317 if reflect:
318 pt = -pt + 2 * cm + 2 * n * np.dot(pt - cm, n)
319 pn.append(pt)
321 n_use.extend(n_use)
322 p_use.extend(pn)
324 # In the case of an uneven number of
325 # atoms we need to add one extra
326 for n in nu:
327 if nu[n] % 2 == 0:
328 continue
329 while sum(n_use == n) > nu[n]:
330 for i in range(int(len(n_use) / 2), len(n_use)):
331 if n_use[i] == n:
332 del p_use[i]
333 del n_use[i]
334 break
335 assert sum(n_use == n) == nu[n]
337 # Make sure we have the correct number of atoms
338 # and rearrange the atoms so they are in the right order
339 for i in range(len(n_use)):
340 if num[i] == n_use[i]:
341 continue
342 for j in range(i + 1, len(n_use)):
343 if n_use[j] == num[i]:
344 tn = n_use[i]
345 tp = p_use[i]
346 n_use[i] = n_use[j]
347 p_use[i] = p_use[j]
348 p_use[j] = tp
349 n_use[j] = tn
351 # Finally we check that nothing is too close in the end product.
352 cand = Atoms(num, p_use, cell=slab.get_cell(), pbc=slab.get_pbc())
354 tc = atoms_too_close(cand, self.blmin)
355 if tc:
356 continue
357 tc = atoms_too_close_two_sets(slab, cand, self.blmin)
359 if not changed and counter > n_tries // 2:
360 reflect = not reflect
361 changed = True
363 tot = slab + cand
365 if counter == n_tries:
366 return None
367 return tot
370class StrainMutation(OffspringCreator):
371 """ Mutates a candidate by applying a randomly generated strain.
373 For more information, see also:
375 * :doi:`Glass, Oganov, Hansen, Comp. Phys. Comm. 175 (2006) 713-720
376 <10.1016/j.cpc.2006.07.020>`
378 * :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
379 <10.1016/j.cpc.2010.07.048>`
381 After initialization of the mutation, a scaling volume
382 (to which each mutated structure is scaled before checking the
383 constraints) is typically generated from the population,
384 which is then also occasionally updated in the course of the
385 GA run.
387 Parameters:
389 blmin: dict
390 The closest allowed interatomic distances on the form:
391 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
393 cellbounds: ase.ga.utilities.CellBounds instance
394 Describes limits on the cell shape, see
395 :class:`~ase.ga.utilities.CellBounds`.
397 stddev: float
398 Standard deviation used in the generation of the
399 strain matrix elements.
401 number_of_variable_cell_vectors: int (default 3)
402 The number of variable cell vectors (1, 2 or 3).
403 To keep things simple, it is the 'first' vectors which
404 will be treated as variable, i.e. the 'a' vector in the
405 univariate case, the 'a' and 'b' vectors in the bivariate
406 case, etc.
408 use_tags: boolean
409 Whether to use the atomic tags to preserve molecular identity.
411 rng: Random number generator
412 By default numpy.random.
413 """
415 def __init__(self, blmin, cellbounds=None, stddev=0.7,
416 number_of_variable_cell_vectors=3, use_tags=False,
417 rng=np.random, verbose=False):
418 OffspringCreator.__init__(self, verbose, rng=rng)
419 self.blmin = blmin
420 self.cellbounds = cellbounds
421 self.stddev = stddev
422 self.number_of_variable_cell_vectors = number_of_variable_cell_vectors
423 self.use_tags = use_tags
425 self.scaling_volume = None
426 self.descriptor = 'StrainMutation'
427 self.min_inputs = 1
429 def update_scaling_volume(self, population, w_adapt=0.5, n_adapt=0):
430 """Function to initialize or update the scaling volume in a GA run.
432 w_adapt: weight of the new vs the old scaling volume
434 n_adapt: number of best candidates in the population that
435 are used to calculate the new scaling volume
436 """
437 if not n_adapt:
438 # if not set, take best 20% of the population
439 n_adapt = int(np.ceil(0.2 * len(population)))
440 v_new = np.mean([a.get_volume() for a in population[:n_adapt]])
442 if not self.scaling_volume:
443 self.scaling_volume = v_new
444 else:
445 volumes = [self.scaling_volume, v_new]
446 weights = [1 - w_adapt, w_adapt]
447 self.scaling_volume = np.average(volumes, weights=weights)
449 def get_new_individual(self, parents):
450 f = parents[0]
452 indi = self.mutate(f)
453 if indi is None:
454 return indi, 'mutation: strain'
456 indi = self.initialize_individual(f, indi)
457 indi.info['data']['parents'] = [f.info['confid']]
459 return self.finalize_individual(indi), 'mutation: strain'
461 def mutate(self, atoms):
462 """ Does the actual mutation. """
463 cell_ref = atoms.get_cell()
464 pos_ref = atoms.get_positions()
466 if self.scaling_volume is None:
467 # The scaling_volume has not been set (yet),
468 # so we give it the same volume as the parent
469 vol_ref = atoms.get_volume()
470 else:
471 vol_ref = self.scaling_volume
473 if self.use_tags:
474 tags = atoms.get_tags()
475 gather_atoms_by_tag(atoms)
476 pos = atoms.get_positions()
478 mutant = atoms.copy()
480 count = 0
481 too_close = True
482 maxcount = 1000
483 while too_close and count < maxcount:
484 count += 1
486 # generating the strain matrix:
487 strain = np.identity(3)
488 for i in range(self.number_of_variable_cell_vectors):
489 for j in range(i + 1):
490 r = self.rng.normal(loc=0., scale=self.stddev)
491 if i == j:
492 strain[i, j] += r
493 else:
494 epsilon = 0.5 * r
495 strain[i, j] += epsilon
496 strain[j, i] += epsilon
498 # applying the strain:
499 cell_new = np.dot(strain, cell_ref)
501 # convert the submatrix with the variable cell vectors
502 # to a lower triangular form
503 cell_new = calc_rotated_cell(cell_new)
504 for i in range(self.number_of_variable_cell_vectors, 3):
505 cell_new[i] = cell_ref[i]
507 cell_new = Cell(cell_new)
509 # volume scaling:
510 if self.number_of_variable_cell_vectors > 0:
511 scaling = vol_ref / cell_new.volume
512 scaling **= 1. / self.number_of_variable_cell_vectors
513 cell_new[:self.number_of_variable_cell_vectors] *= scaling
515 # check cell dimensions:
516 if not self.cellbounds.is_within_bounds(cell_new):
517 continue
519 # ensure non-variable cell vectors are indeed unchanged
520 for i in range(self.number_of_variable_cell_vectors, 3):
521 assert np.allclose(cell_new[i], cell_ref[i])
523 # check that the volume is correct
524 assert np.isclose(vol_ref, cell_new.volume)
526 # apply the new unit cell and scale
527 # the atomic positions accordingly
528 mutant.set_cell(cell_ref, scale_atoms=False)
530 if self.use_tags:
531 transfo = np.linalg.solve(cell_ref, cell_new)
532 for tag in np.unique(tags):
533 select = np.where(tags == tag)
534 cop = np.mean(pos[select], axis=0)
535 disp = np.dot(cop, transfo) - cop
536 mutant.positions[select] += disp
537 else:
538 mutant.set_positions(pos_ref)
540 mutant.set_cell(cell_new, scale_atoms=not self.use_tags)
541 mutant.wrap()
543 # check the interatomic distances
544 too_close = atoms_too_close(mutant, self.blmin,
545 use_tags=self.use_tags)
547 if count == maxcount:
548 mutant = None
550 return mutant
553class PermuStrainMutation(CombinationMutation):
554 """Combination of PermutationMutation and StrainMutation.
556 For more information, see also:
558 * :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
559 <10.1016/j.cpc.2010.07.048>`
561 Parameters:
563 permutationmutation: OffspringCreator instance
564 A mutation that permutes atom types.
566 strainmutation: OffspringCreator instance
567 A mutation that mutates by straining.
568 """
570 def __init__(self, permutationmutation, strainmutation, verbose=False):
571 super().__init__(permutationmutation,
572 strainmutation,
573 verbose=verbose)
574 self.descriptor = 'permustrain'
577class RotationalMutation(OffspringCreator):
578 """Mutates a candidate by applying random rotations
579 to multi-atom moieties in the structure (atoms with
580 the same tag are considered part of one such moiety).
582 Only performs whole-molecule rotations, no internal
583 rotations.
585 For more information, see also:
587 * `Zhu Q., Oganov A.R., Glass C.W., Stokes H.T,
588 Acta Cryst. (2012), B68, 215-226.`__
590 __ https://dx.doi.org/10.1107/S0108768112017466
592 Parameters:
594 blmin: dict
595 The closest allowed interatomic distances on the form:
596 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
598 n_top: int or None
599 The number of atoms to optimize (None = include all).
601 fraction: float
602 Fraction of the moieties to be rotated.
604 tags: None or list of integers
605 Specifies, respectively, whether all moieties or only those
606 with matching tags are eligible for rotation.
608 min_angle: float
609 Minimal angle (in radians) for each rotation;
610 should lie in the interval [0, pi].
612 test_dist_to_slab: boolean
613 Whether also the distances to the slab
614 should be checked to satisfy the blmin.
616 rng: Random number generator
617 By default numpy.random.
618 """
620 def __init__(self, blmin, n_top=None, fraction=0.33, tags=None,
621 min_angle=1.57, test_dist_to_slab=True, rng=np.random,
622 verbose=False):
623 OffspringCreator.__init__(self, verbose, rng=rng)
624 self.blmin = blmin
625 self.n_top = n_top
626 self.fraction = fraction
627 self.tags = tags
628 self.min_angle = min_angle
629 self.test_dist_to_slab = test_dist_to_slab
630 self.descriptor = 'RotationalMutation'
631 self.min_inputs = 1
633 def get_new_individual(self, parents):
634 f = parents[0]
636 indi = self.mutate(f)
637 if indi is None:
638 return indi, 'mutation: rotational'
640 indi = self.initialize_individual(f, indi)
641 indi.info['data']['parents'] = [f.info['confid']]
643 return self.finalize_individual(indi), 'mutation: rotational'
645 def mutate(self, atoms):
646 """Does the actual mutation."""
647 N = len(atoms) if self.n_top is None else self.n_top
648 slab = atoms[:len(atoms) - N]
649 atoms = atoms[-N:]
651 mutant = atoms.copy()
652 gather_atoms_by_tag(mutant)
653 pos = mutant.get_positions()
654 tags = mutant.get_tags()
655 eligible_tags = tags if self.tags is None else self.tags
657 indices = {}
658 for tag in np.unique(tags):
659 hits = np.where(tags == tag)[0]
660 if len(hits) > 1 and tag in eligible_tags:
661 indices[tag] = hits
663 n_rot = int(np.ceil(len(indices) * self.fraction))
664 chosen_tags = self.rng.choice(list(indices.keys()), size=n_rot,
665 replace=False)
667 too_close = True
668 count = 0
669 maxcount = 10000
670 while too_close and count < maxcount:
671 newpos = np.copy(pos)
672 for tag in chosen_tags:
673 p = np.copy(newpos[indices[tag]])
674 cop = np.mean(p, axis=0)
676 if len(p) == 2:
677 line = (p[1] - p[0]) / np.linalg.norm(p[1] - p[0])
678 while True:
679 axis = self.rng.random(3)
680 axis /= np.linalg.norm(axis)
681 a = np.arccos(np.dot(axis, line))
682 if np.pi / 4 < a < np.pi * 3 / 4:
683 break
684 else:
685 axis = self.rng.random(3)
686 axis /= np.linalg.norm(axis)
688 angle = self.min_angle
689 angle += 2 * (np.pi - self.min_angle) * self.rng.random()
691 m = get_rotation_matrix(axis, angle)
692 newpos[indices[tag]] = np.dot(m, (p - cop).T).T + cop
694 mutant.set_positions(newpos)
695 mutant.wrap()
696 too_close = atoms_too_close(mutant, self.blmin, use_tags=True)
697 count += 1
699 if not too_close and self.test_dist_to_slab:
700 too_close = atoms_too_close_two_sets(slab, mutant, self.blmin)
702 if count == maxcount:
703 mutant = None
704 else:
705 mutant = slab + mutant
707 return mutant
710class RattleRotationalMutation(CombinationMutation):
711 """Combination of RattleMutation and RotationalMutation.
713 Parameters:
715 rattlemutation: OffspringCreator instance
716 A mutation that rattles atoms.
718 rotationalmutation: OffspringCreator instance
719 A mutation that rotates moieties.
720 """
722 def __init__(self, rattlemutation, rotationalmutation, verbose=False):
723 super().__init__(rattlemutation,
724 rotationalmutation,
725 verbose=verbose)
726 self.descriptor = 'rattlerotational'