Coverage for /builds/kinetik161/ase/ase/ga/particle_crossovers.py: 73.91%
115 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1from itertools import chain
3import numpy as np
5from ase import Atoms
6from ase.ga.offspring_creator import OffspringCreator
9class Crossover(OffspringCreator):
10 """Base class for all particle crossovers.
12 Originally intended for medium sized particles
14 Do not call this class directly."""
16 def __init__(self, rng=np.random):
17 OffspringCreator.__init__(self, rng=rng)
18 self.descriptor = 'Crossover'
19 self.min_inputs = 2
22class CutSpliceCrossover(Crossover):
23 """Crossover that cuts two particles through a plane in space and
24 merges two halfes from different particles together.
26 Implementation of the method presented in:
27 D. M. Deaven and K. M. Ho, Phys. Rev. Lett., 75, 2, 288-291 (1995)
29 It keeps the correct composition by randomly assigning elements in
30 the new particle. If some of the atoms in the two particle halves
31 are too close, the halves are moved away from each other perpendicular
32 to the cutting plane.
34 Parameters:
36 blmin: dictionary of minimum distance between atomic numbers.
37 e.g. {(28,29): 1.5}
39 keep_composition: boolean that signifies if the composition should
40 be the same as in the parents.
42 rng: Random number generator
43 By default numpy.random.
44 """
46 def __init__(self, blmin, keep_composition=True, rng=np.random):
47 Crossover.__init__(self, rng=rng)
48 self.blmin = blmin
49 self.keep_composition = keep_composition
50 self.descriptor = 'CutSpliceCrossover'
52 def get_new_individual(self, parents):
53 f, m = parents
55 indi = self.initialize_individual(f)
56 indi.info['data']['parents'] = [i.info['confid'] for i in parents]
58 theta = self.rng.random() * 2 * np.pi # 0,2pi
59 phi = self.rng.random() * np.pi # 0,pi
60 e = np.array((np.sin(phi) * np.cos(theta),
61 np.sin(theta) * np.sin(phi),
62 np.cos(phi)))
63 eps = 0.0001
65 f.translate(-f.get_center_of_mass())
66 m.translate(-m.get_center_of_mass())
68 # Get the signed distance to the cutting plane
69 # We want one side from f and the other side from m
70 fmap = [np.dot(x, e) for x in f.get_positions()]
71 mmap = [-np.dot(x, e) for x in m.get_positions()]
72 ain = sorted([i for i in chain(fmap, mmap) if i > 0],
73 reverse=True)
74 aout = sorted([i for i in chain(fmap, mmap) if i < 0],
75 reverse=True)
77 off = len(ain) - len(f)
79 # Translating f and m to get the correct number of atoms
80 # in the offspring
81 if off < 0:
82 # too few
83 # move f and m away from the plane
84 dist = (abs(aout[abs(off) - 1]) + abs(aout[abs(off)])) * .5
85 f.translate(e * dist)
86 m.translate(-e * dist)
87 elif off > 0:
88 # too many
89 # move f and m towards the plane
90 dist = (abs(ain[-off - 1]) + abs(ain[-off])) * .5
91 f.translate(-e * dist)
92 m.translate(e * dist)
93 if off != 0 and dist == 0:
94 # Exactly same position => we continue with the wrong number
95 # of atoms. What should be done? Fail or return None or
96 # remove one of the two atoms with exactly the same position.
97 pass
99 # Determine the contributing parts from f and m
100 tmpf, tmpm = Atoms(), Atoms()
101 for atom in f:
102 if np.dot(atom.position, e) > 0:
103 atom.tag = 1
104 tmpf.append(atom)
105 for atom in m:
106 if np.dot(atom.position, e) < 0:
107 atom.tag = 2
108 tmpm.append(atom)
110 # Check that the correct composition is employed
111 if self.keep_composition:
112 opt_sm = sorted(f.numbers)
113 tmpf_numbers = list(tmpf.numbers)
114 tmpm_numbers = list(tmpm.numbers)
115 cur_sm = sorted(tmpf_numbers + tmpm_numbers)
116 # correct_by: dictionary that specifies how many
117 # of the atom_numbers should be removed (a negative number)
118 # or added (a positive number)
119 correct_by = {j: opt_sm.count(j) for j in set(opt_sm)}
120 for n in cur_sm:
121 correct_by[n] -= 1
122 correct_in = tmpf if self.rng.choice([0, 1]) else tmpm
123 to_add, to_rem = [], []
124 for num, amount in correct_by.items():
125 if amount > 0:
126 to_add.extend([num] * amount)
127 elif amount < 0:
128 to_rem.extend([num] * abs(amount))
129 for add, rem in zip(to_add, to_rem):
130 tbc = [a.index for a in correct_in if a.number == rem]
131 if len(tbc) == 0:
132 pass
133 ai = self.rng.choice(tbc)
134 correct_in[ai].number = add
136 # Move the contributing apart if any distance is below blmin
137 maxl = 0.
138 for sv, min_dist in self.get_vectors_below_min_dist(tmpf + tmpm):
139 lsv = np.linalg.norm(sv) # length of shortest vector
140 d = [-np.dot(e, sv)] * 2
141 d[0] += np.sqrt(np.dot(e, sv)**2 - lsv**2 + min_dist**2)
142 d[1] -= np.sqrt(np.dot(e, sv)**2 - lsv**2 + min_dist**2)
143 L = sorted([abs(i) for i in d])[0] / 2. + eps
144 if L > maxl:
145 maxl = L
146 tmpf.translate(e * maxl)
147 tmpm.translate(-e * maxl)
149 # Put the two parts together
150 for atom in chain(tmpf, tmpm):
151 indi.append(atom)
153 parent_message = ':Parents {} {}'.format(f.info['confid'],
154 m.info['confid'])
155 return (self.finalize_individual(indi),
156 self.descriptor + parent_message)
158 def get_numbers(self, atoms):
159 """Returns the atomic numbers of the atoms object using only
160 the elements defined in self.elements"""
161 ac = atoms.copy()
162 if self.elements is not None:
163 del ac[[a.index for a in ac
164 if a.symbol in self.elements]]
165 return ac.numbers
167 def get_vectors_below_min_dist(self, atoms):
168 """Generator function that returns each vector (between atoms)
169 that is shorter than the minimum distance for those atom types
170 (set during the initialization in blmin)."""
171 norm = np.linalg.norm
172 ap = atoms.get_positions()
173 an = atoms.numbers
174 for i in range(len(atoms)):
175 pos = atoms[i].position
176 for j, d in enumerate([norm(k - pos) for k in ap[i:]]):
177 if d == 0:
178 continue
179 min_dist = self.blmin[tuple(sorted((an[i], an[j + i])))]
180 if d < min_dist:
181 yield atoms[i].position - atoms[j + i].position, min_dist
183 def get_shortest_dist_vector(self, atoms):
184 norm = np.linalg.norm
185 mind = 10000.
186 ap = atoms.get_positions()
187 for i in range(len(atoms)):
188 pos = atoms[i].position
189 for j, d in enumerate([norm(k - pos) for k in ap[i:]]):
190 if d == 0:
191 continue
192 if d < mind:
193 mind = d
194 lowpair = (i, j + i)
195 return atoms[lowpair[0]].position - atoms[lowpair[1]].position