Coverage for /builds/debichem-team/python-ase/ase/ga/standardmutations.py : 73.32%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""A collection of mutations that can be used."""
2import numpy as np
3from math import cos, sin, pi
4from ase.calculators.lammpslib import convert_cell
5from ase.ga.utilities import (atoms_too_close,
6 atoms_too_close_two_sets,
7 gather_atoms_by_tag,
8 get_rotation_matrix)
9from ase.ga.offspring_creator import OffspringCreator, CombinationMutation
10from ase import Atoms
13class RattleMutation(OffspringCreator):
14 """An implementation of the rattle mutation as described in:
16 R.L. Johnston Dalton Transactions, Vol. 22,
17 No. 22. (2003), pp. 4193-4207
19 Parameters:
21 blmin: Dictionary defining the minimum distance between atoms
22 after the rattle.
24 n_top: Number of atoms optimized by the GA.
26 rattle_strength: Strength with which the atoms are moved.
28 rattle_prop: The probability with which each atom is rattled.
30 test_dist_to_slab: whether to also make sure that the distances
31 between the atoms and the slab satisfy the blmin.
33 use_tags: if True, the atomic tags will be used to preserve
34 molecular identity. Same-tag atoms will then be
35 displaced collectively, so that the internal
36 geometry is preserved.
38 rng: Random number generator
39 By default numpy.random.
40 """
41 def __init__(self, blmin, n_top, rattle_strength=0.8,
42 rattle_prop=0.4, test_dist_to_slab=True, use_tags=False,
43 verbose=False, rng=np.random):
44 OffspringCreator.__init__(self, verbose, rng=rng)
45 self.blmin = blmin
46 self.n_top = n_top
47 self.rattle_strength = rattle_strength
48 self.rattle_prop = rattle_prop
49 self.test_dist_to_slab = test_dist_to_slab
50 self.use_tags = use_tags
52 self.descriptor = 'RattleMutation'
53 self.min_inputs = 1
55 def get_new_individual(self, parents):
56 f = parents[0]
58 indi = self.mutate(f)
59 if indi is None:
60 return indi, 'mutation: rattle'
62 indi = self.initialize_individual(f, indi)
63 indi.info['data']['parents'] = [f.info['confid']]
65 return self.finalize_individual(indi), 'mutation: rattle'
67 def mutate(self, atoms):
68 """Does the actual mutation."""
69 N = len(atoms) if self.n_top is None else self.n_top
70 slab = atoms[:len(atoms) - N]
71 atoms = atoms[-N:]
72 tags = atoms.get_tags() if self.use_tags else np.arange(N)
73 pos_ref = atoms.get_positions()
74 num = atoms.get_atomic_numbers()
75 cell = atoms.get_cell()
76 pbc = atoms.get_pbc()
77 st = 2. * self.rattle_strength
79 count = 0
80 maxcount = 1000
81 too_close = True
82 while too_close and count < maxcount:
83 count += 1
84 pos = pos_ref.copy()
85 ok = False
86 for tag in np.unique(tags):
87 select = np.where(tags == tag)
88 if self.rng.rand() < self.rattle_prop:
89 ok = True
90 r = self.rng.rand(3)
91 pos[select] += st * (r - 0.5)
93 if not ok:
94 # Nothing got rattled
95 continue
97 top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
98 too_close = atoms_too_close(
99 top, self.blmin, use_tags=self.use_tags)
100 if not too_close and self.test_dist_to_slab:
101 too_close = atoms_too_close_two_sets(top, slab, self.blmin)
103 if count == maxcount:
104 return None
106 mutant = slab + top
107 return mutant
110class PermutationMutation(OffspringCreator):
111 """Mutation that permutes a percentage of the atom types in the cluster.
113 Parameters:
115 n_top: Number of atoms optimized by the GA.
117 probability: The probability with which an atom is permuted.
119 test_dist_to_slab: whether to also make sure that the distances
120 between the atoms and the slab satisfy the blmin.
122 use_tags: if True, the atomic tags will be used to preserve
123 molecular identity. Permutations will then happen
124 at the molecular level, i.e. swapping the center-of-
125 positions of two moieties while preserving their
126 internal geometries.
128 blmin: Dictionary defining the minimum distance between atoms
129 after the permutation. If equal to None (the default),
130 no such check is performed.
132 rng: Random number generator
133 By default numpy.random.
134 """
135 def __init__(self, n_top, probability=0.33, test_dist_to_slab=True,
136 use_tags=False, blmin=None, rng=np.random, verbose=False):
137 OffspringCreator.__init__(self, verbose, rng=rng)
138 self.n_top = n_top
139 self.probability = probability
140 self.test_dist_to_slab = test_dist_to_slab
141 self.use_tags = use_tags
142 self.blmin = blmin
144 self.descriptor = 'PermutationMutation'
145 self.min_inputs = 1
147 def get_new_individual(self, parents):
148 f = parents[0]
150 indi = self.mutate(f)
151 if indi is None:
152 return indi, 'mutation: permutation'
154 indi = self.initialize_individual(f, indi)
155 indi.info['data']['parents'] = [f.info['confid']]
157 return self.finalize_individual(indi), 'mutation: permutation'
159 def mutate(self, atoms):
160 """Does the actual mutation."""
161 N = len(atoms) if self.n_top is None else self.n_top
162 slab = atoms[:len(atoms) - N]
163 atoms = atoms[-N:]
164 if self.use_tags:
165 gather_atoms_by_tag(atoms)
166 tags = atoms.get_tags() if self.use_tags else np.arange(N)
167 pos_ref = atoms.get_positions()
168 num = atoms.get_atomic_numbers()
169 cell = atoms.get_cell()
170 pbc = atoms.get_pbc()
171 symbols = atoms.get_chemical_symbols()
173 unique_tags = np.unique(tags)
174 n = len(unique_tags)
175 swaps = int(np.ceil(n * self.probability / 2.))
177 sym = []
178 for tag in unique_tags:
179 indices = np.where(tags == tag)[0]
180 s = ''.join([symbols[j] for j in indices])
181 sym.append(s)
183 assert len(np.unique(sym)) > 1, \
184 'Permutations with one atom (or molecule) type is not valid'
186 count = 0
187 maxcount = 1000
188 too_close = True
189 while too_close and count < maxcount:
190 count += 1
191 pos = pos_ref.copy()
192 for _ in range(swaps):
193 i = j = 0
194 while sym[i] == sym[j]:
195 i = self.rng.randint(0, high=n)
196 j = self.rng.randint(0, high=n)
197 ind1 = np.where(tags == i)
198 ind2 = np.where(tags == j)
199 cop1 = np.mean(pos[ind1], axis=0)
200 cop2 = np.mean(pos[ind2], axis=0)
201 pos[ind1] += cop2 - cop1
202 pos[ind2] += cop1 - cop2
204 top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
205 if self.blmin is None:
206 too_close = False
207 else:
208 too_close = atoms_too_close(
209 top, self.blmin, use_tags=self.use_tags)
210 if not too_close and self.test_dist_to_slab:
211 too_close = atoms_too_close_two_sets(top, slab, self.blmin)
213 if count == maxcount:
214 return None
216 mutant = slab + top
217 return mutant
220class MirrorMutation(OffspringCreator):
221 """A mirror mutation, as described in
222 TO BE PUBLISHED.
224 This mutation mirrors half of the cluster in a
225 randomly oriented cutting plane discarding the other half.
227 Parameters:
229 blmin: Dictionary defining the minimum allowed
230 distance between atoms.
232 n_top: Number of atoms the GA optimizes.
234 reflect: Defines if the mirrored half is also reflected
235 perpendicular to the mirroring plane.
237 rng: Random number generator
238 By default numpy.random.
239 """
240 def __init__(self, blmin, n_top, reflect=False, rng=np.random,
241 verbose=False):
242 OffspringCreator.__init__(self, verbose, rng=rng)
243 self.blmin = blmin
244 self.n_top = n_top
245 self.reflect = reflect
247 self.descriptor = 'MirrorMutation'
248 self.min_inputs = 1
250 def get_new_individual(self, parents):
251 f = parents[0]
253 indi = self.mutate(f)
254 if indi is None:
255 return indi, 'mutation: mirror'
257 indi = self.initialize_individual(f, indi)
258 indi.info['data']['parents'] = [f.info['confid']]
260 return self.finalize_individual(indi), 'mutation: mirror'
262 def mutate(self, atoms):
263 """ Do the mutation of the atoms input. """
264 reflect = self.reflect
265 tc = True
266 slab = atoms[0:len(atoms) - self.n_top]
267 top = atoms[len(atoms) - self.n_top: len(atoms)]
268 num = top.numbers
269 unique_types = list(set(num))
270 nu = dict()
271 for u in unique_types:
272 nu[u] = sum(num == u)
274 n_tries = 1000
275 counter = 0
276 changed = False
278 while tc and counter < n_tries:
279 counter += 1
280 cand = top.copy()
281 pos = cand.get_positions()
283 cm = np.average(top.get_positions(), axis=0)
285 # first select a randomly oriented cutting plane
286 theta = pi * self.rng.rand()
287 phi = 2. * pi * self.rng.rand()
288 n = (cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta))
289 n = np.array(n)
291 # Calculate all atoms signed distance to the cutting plane
292 D = []
293 for (i, p) in enumerate(pos):
294 d = np.dot(p - cm, n)
295 D.append((i, d))
297 # Sort the atoms by their signed distance
298 D.sort(key=lambda x: x[1])
299 nu_taken = dict()
301 # Select half of the atoms needed for a full cluster
302 p_use = []
303 n_use = []
304 for (i, d) in D:
305 if num[i] not in nu_taken.keys():
306 nu_taken[num[i]] = 0
307 if nu_taken[num[i]] < nu[num[i]] / 2.:
308 p_use.append(pos[i])
309 n_use.append(num[i])
310 nu_taken[num[i]] += 1
312 # calculate the mirrored position and add these.
313 pn = []
314 for p in p_use:
315 pt = p - 2. * np.dot(p - cm, n) * n
316 if reflect:
317 pt = -pt + 2 * cm + 2 * n * np.dot(pt - cm, n)
318 pn.append(pt)
320 n_use.extend(n_use)
321 p_use.extend(pn)
323 # In the case of an uneven number of
324 # atoms we need to add one extra
325 for n in nu.keys():
326 if nu[n] % 2 == 0:
327 continue
328 while sum(n_use == n) > nu[n]:
329 for i in range(int(len(n_use) / 2), len(n_use)):
330 if n_use[i] == n:
331 del p_use[i]
332 del n_use[i]
333 break
334 assert sum(n_use == n) == nu[n]
336 # Make sure we have the correct number of atoms
337 # and rearrange the atoms so they are in the right order
338 for i in range(len(n_use)):
339 if num[i] == n_use[i]:
340 continue
341 for j in range(i + 1, len(n_use)):
342 if n_use[j] == num[i]:
343 tn = n_use[i]
344 tp = p_use[i]
345 n_use[i] = n_use[j]
346 p_use[i] = p_use[j]
347 p_use[j] = tp
348 n_use[j] = tn
350 # Finally we check that nothing is too close in the end product.
351 cand = Atoms(num, p_use, cell=slab.get_cell(), pbc=slab.get_pbc())
353 tc = atoms_too_close(cand, self.blmin)
354 if tc:
355 continue
356 tc = atoms_too_close_two_sets(slab, cand, self.blmin)
358 if not changed and counter > n_tries // 2:
359 reflect = not reflect
360 changed = True
362 tot = slab + cand
364 if counter == n_tries:
365 return None
366 return tot
369class StrainMutation(OffspringCreator):
370 """ Mutates a candidate by applying a randomly generated strain.
372 For more information, see also:
374 * `Glass, Oganov, Hansen, Comp. Phys. Comm. 175 (2006) 713-720`__
376 __ https://doi.org/10.1016/j.cpc.2006.07.020
378 * `Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387`__
380 __ https://doi.org/10.1016/j.cpc.2010.07.048
382 After initialization of the mutation, a scaling volume
383 (to which each mutated structure is scaled before checking the
384 constraints) is typically generated from the population,
385 which is then also occasionally updated in the course of the
386 GA run.
388 Parameters:
390 blmin: dict
391 The closest allowed interatomic distances on the form:
392 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
394 cellbounds: ase.ga.utilities.CellBounds instance
395 Describes limits on the cell shape, see
396 :class:`~ase.ga.utilities.CellBounds`.
398 stddev: float
399 Standard deviation used in the generation of the
400 strain matrix elements.
402 number_of_variable_cell_vectors: int (default 3)
403 The number of variable cell vectors (1, 2 or 3).
404 To keep things simple, it is the 'first' vectors which
405 will be treated as variable, i.e. the 'a' vector in the
406 univariate case, the 'a' and 'b' vectors in the bivariate
407 case, etc.
409 use_tags: boolean
410 Whether to use the atomic tags to preserve molecular identity.
412 rng: Random number generator
413 By default numpy.random.
414 """
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()
465 vol_ref = atoms.get_volume()
467 if self.use_tags:
468 tags = atoms.get_tags()
469 gather_atoms_by_tag(atoms)
470 pos = atoms.get_positions()
472 mutant = atoms.copy()
474 count = 0
475 too_close = True
476 maxcount = 1000
477 while too_close and count < maxcount:
478 count += 1
480 # generating the strain matrix:
481 strain = np.identity(3)
482 for i in range(self.number_of_variable_cell_vectors):
483 for j in range(i + 1):
484 r = self.rng.normal(loc=0., scale=self.stddev)
485 if i == j:
486 strain[i, j] += r
487 else:
488 epsilon = 0.5 * r
489 strain[i, j] += epsilon
490 strain[j, i] += epsilon
492 # applying the strain:
493 cell_new = np.dot(strain, cell_ref)
495 # convert to lower triangular form
496 cell_new = convert_cell(cell_new)[0].T
498 # volume scaling:
499 if self.number_of_variable_cell_vectors > 0:
500 volume = abs(np.linalg.det(cell_new))
501 if self.scaling_volume is None:
502 # The scaling_volume has not been set (yet),
503 # so we give it the same volume as the parent
504 scaling = vol_ref / volume
505 else:
506 scaling = self.scaling_volume / volume
507 scaling **= 1. / self.number_of_variable_cell_vectors
508 cell_new[:self.number_of_variable_cell_vectors] *= scaling
510 # check cell dimensions:
511 if not self.cellbounds.is_within_bounds(cell_new):
512 continue
514 # ensure non-variable cell vectors are indeed unchanged
515 for i in range(self.number_of_variable_cell_vectors, 3):
516 assert np.allclose(cell_new[i], cell_ref[i])
518 # apply the new unit cell and scale
519 # the atomic positions accordingly
520 mutant.set_cell(cell_ref, scale_atoms=False)
522 if self.use_tags:
523 transfo = np.linalg.solve(cell_ref, cell_new)
524 for tag in np.unique(tags):
525 select = np.where(tags == tag)
526 cop = np.mean(pos[select], axis=0)
527 disp = np.dot(cop, transfo) - cop
528 mutant.positions[select] += disp
529 else:
530 mutant.set_positions(pos_ref)
532 mutant.set_cell(cell_new, scale_atoms=not self.use_tags)
533 mutant.wrap()
535 # check the interatomic distances
536 too_close = atoms_too_close(mutant, self.blmin,
537 use_tags=self.use_tags)
539 if count == maxcount:
540 mutant = None
542 return mutant
545class PermuStrainMutation(CombinationMutation):
546 """Combination of PermutationMutation and StrainMutation.
548 For more information, see also:
550 * `Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387`__
552 __ https://doi.org/10.1016/j.cpc.2010.07.048
554 Parameters:
556 permutationmutation: OffspringCreator instance
557 A mutation that permutes atom types.
559 strainmutation: OffspringCreator instance
560 A mutation that mutates by straining.
561 """
562 def __init__(self, permutationmutation, strainmutation, verbose=False):
563 super(PermuStrainMutation, self).__init__(permutationmutation,
564 strainmutation,
565 verbose=verbose)
566 self.descriptor = 'permustrain'
569class RotationalMutation(OffspringCreator):
570 """Mutates a candidate by applying random rotations
571 to multi-atom moieties in the structure (atoms with
572 the same tag are considered part of one such moiety).
574 Only performs whole-molecule rotations, no internal
575 rotations.
577 For more information, see also:
579 * `Zhu Q., Oganov A.R., Glass C.W., Stokes H.T,
580 Acta Cryst. (2012), B68, 215-226.`__
582 __ https://dx.doi.org/10.1107/S0108768112017466
584 Parameters:
586 blmin: dict
587 The closest allowed interatomic distances on the form:
588 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers.
590 n_top: int or None
591 The number of atoms to optimize (None = include all).
593 fraction: float
594 Fraction of the moieties to be rotated.
596 tags: None or list of integers
597 Specifies, respectively, whether all moieties or only those
598 with matching tags are eligible for rotation.
600 min_angle: float
601 Minimal angle (in radians) for each rotation;
602 should lie in the interval [0, pi].
604 test_dist_to_slab: boolean
605 Whether also the distances to the slab
606 should be checked to satisfy the blmin.
608 rng: Random number generator
609 By default numpy.random.
610 """
611 def __init__(self, blmin, n_top=None, fraction=0.33, tags=None,
612 min_angle=1.57, test_dist_to_slab=True, rng=np.random,
613 verbose=False):
614 OffspringCreator.__init__(self, verbose, rng=rng)
615 self.blmin = blmin
616 self.n_top = n_top
617 self.fraction = fraction
618 self.tags = tags
619 self.min_angle = min_angle
620 self.test_dist_to_slab = test_dist_to_slab
621 self.descriptor = 'RotationalMutation'
622 self.min_inputs = 1
624 def get_new_individual(self, parents):
625 f = parents[0]
627 indi = self.mutate(f)
628 if indi is None:
629 return indi, 'mutation: rotational'
631 indi = self.initialize_individual(f, indi)
632 indi.info['data']['parents'] = [f.info['confid']]
634 return self.finalize_individual(indi), 'mutation: rotational'
636 def mutate(self, atoms):
637 """Does the actual mutation."""
638 N = len(atoms) if self.n_top is None else self.n_top
639 slab = atoms[:len(atoms) - N]
640 atoms = atoms[-N:]
642 mutant = atoms.copy()
643 gather_atoms_by_tag(mutant)
644 pos = mutant.get_positions()
645 tags = mutant.get_tags()
646 eligible_tags = tags if self.tags is None else self.tags
648 indices = {}
649 for tag in np.unique(tags):
650 hits = np.where(tags == tag)[0]
651 if len(hits) > 1 and tag in eligible_tags:
652 indices[tag] = hits
654 n_rot = int(np.ceil(len(indices) * self.fraction))
655 chosen_tags = self.rng.choice(list(indices.keys()), size=n_rot,
656 replace=False)
658 too_close = True
659 count = 0
660 maxcount = 10000
661 while too_close and count < maxcount:
662 newpos = np.copy(pos)
663 for tag in chosen_tags:
664 p = np.copy(newpos[indices[tag]])
665 cop = np.mean(p, axis=0)
667 if len(p) == 2:
668 line = (p[1] - p[0]) / np.linalg.norm(p[1] - p[0])
669 while True:
670 axis = self.rng.rand(3)
671 axis /= np.linalg.norm(axis)
672 a = np.arccos(np.dot(axis, line))
673 if np.pi / 4 < a < np.pi * 3 / 4:
674 break
675 else:
676 axis = self.rng.rand(3)
677 axis /= np.linalg.norm(axis)
679 angle = self.min_angle
680 angle += 2 * (np.pi - self.min_angle) * self.rng.rand()
682 m = get_rotation_matrix(axis, angle)
683 newpos[indices[tag]] = np.dot(m, (p - cop).T).T + cop
685 mutant.set_positions(newpos)
686 mutant.wrap()
687 too_close = atoms_too_close(mutant, self.blmin, use_tags=True)
688 count += 1
690 if not too_close and self.test_dist_to_slab:
691 too_close = atoms_too_close_two_sets(slab, mutant, self.blmin)
693 if count == maxcount:
694 mutant = None
695 else:
696 mutant = slab + mutant
698 return mutant
701class RattleRotationalMutation(CombinationMutation):
702 """Combination of RattleMutation and RotationalMutation.
704 Parameters:
706 rattlemutation: OffspringCreator instance
707 A mutation that rattles atoms.
709 rotationalmutation: OffspringCreator instance
710 A mutation that rotates moieties.
711 """
712 def __init__(self, rattlemutation, rotationalmutation, verbose=False):
713 super(RattleRotationalMutation, self).__init__(rattlemutation,
714 rotationalmutation,
715 verbose=verbose)
716 self.descriptor = 'rattlerotational'