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