standardmutations_mkb.py 17.3 KB
Newer Older
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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
""" A collection of mutations that can be used. """

import numpy as np
from math import ceil, cos, sin, pi
from ase.ga.utilities import atoms_too_close
from ase.ga.utilities import atoms_too_close_two_sets
from mutations.offspring_creator import OffspringCreator
from ase import Atoms, Atom

def pos_add_sphere(rattle_strength):
    # Rattle within a sphere
    r = rattle_strength * np.random.rand()**(1/3)
    theta = np.random.uniform(low=0, high=2*np.pi)
    phi = np.random.uniform(low=0, high=np.pi)
    pos_add = r * np.array([np.cos(theta)*np.sin(phi),
                            np.sin(theta)*np.sin(phi),
                            np.cos(phi)])
    return pos_add

def pos_add_cylinder(rattle_strength, direction, atoms):
    direction = int(direction)
    cylinder_length = atoms.cell[direction][direction]
    # Rattle within a cylinder
    r = rattle_strength * np.random.rand()**(1/2)
    theta = np.random.uniform(low=0, high=2*np.pi)

    pos_add_axis = np.random.uniform(low=-0.5*cylinder_length, high=0.5*cylinder_length)
    pos_add_radial = r * np.array([np.cos(theta),
                                   np.sin(theta)])
    if direction == 0:
        pos_add = np.array([pos_add_axis,
                            pos_add_radial[0],
                            pos_add_radial[1]])
    elif direction == 1:
        pos_add = r * np.array([pos_add_radial[0],
                                pos_add_axis,
                                pos_add_radial[1]])
    elif direction == 2:
        pos_add = r * np.array([pos_add_radial[0],
                                pos_add_radial[1],
                                pos_add_axis])
    else:
        raise RuntimeError('cylinder direction can only be 0 or 1 or 2')
    return pos_add


class RattleMutation(OffspringCreator):
    """ 
    An implementation of the rattle mutation as described in
    R.L. Johnston Dalton Transactions, Vol. 22,
    No. 22. (2003), pp. 4193-4207

    Parameters:

    blmin: Dictionary defining the minimum distance between atoms
    after the rattle.
    n_top: Number of atoms optimized by the GA.
    rattle_strength: Strength with which the atoms are moved.
    rattle_prop: The probability with which each atom is rattled.
    min_z: minimum allowed z-coordinate for rattled atoms.
    rattle_mode: 'sphere' or 'cylinder'
    """
    def __init__(self, blmin, n_top, rattle_strength=0.8,
                 rattle_prop=0.4, min_z=None, use_max_dist=True, verbose=False,
                 rattle_mode='sphere', cluster=False, cylinder_direction=0,
                 descriptor='RattleMutation'):
        OffspringCreator.__init__(self, verbose)
        self.blmin = blmin
        self.n_top = n_top
        self.rattle_strength = rattle_strength
        self.rattle_prop = rattle_prop
        self.min_z = min_z
        self.use_max_dist = use_max_dist
        self.rattle_mode = rattle_mode
        self.cluster = cluster
        self.cylinder_direction = cylinder_direction
        self.descriptor = descriptor
        self.min_inputs = 1

    def get_new_individual(self, parents, **kwargs):
        f = parents[0]

        indi = self.mutate(f)
        if indi is None:
            return indi, 'mutation: rattle'
            
        indi = self.initialize_individual(f, indi)
        #indi.info['data']['parents'] = [f.info['confid']]  # remove for GP-search

        return self.finalize_individual(indi), 'mutation: rattle'

    def mutate(self, atoms):
        num = atoms.numbers
        atoms_mutated = atoms.copy()
        Natoms = len(atoms)
        Nslab = Natoms - self.n_top
        for i in range(Nslab, Natoms):
            if np.random.random() < self.rattle_prop:
                posi_0 = np.copy(atoms_mutated.positions[i])
                tc = False
                tf = False
                for k in range(100):
                    atoms_mutated.positions[i] = posi_0

                    # Make random displacement of atom i
                    if self.rattle_mode == 'sphere':
                        pos_add = pos_add_sphere(self.rattle_strength)
                    elif self.rattle_mode == 'cylinder':
                        pos_add = pos_add_cylinder(self.rattle_strength, self.cylinder_direction, atoms)
                    atoms_mutated.positions[i] += pos_add

                    # Too low - made to avoid that atom is placed below slab
                    if self.min_z is not None:
                        tl = atoms_mutated.positions[i][2] < self.min_z
                        if tl:
                            continue
                    else:
                        tl = False
                    
                    # too close
                    tc = atoms_too_close_two_sets(atoms_mutated[:i] + atoms_mutated[i+1:],
                                         atoms_mutated[i], self.blmin)
                    # too far
                    if self.use_max_dist:
                        if self.cluster:
                            index_other = np.delete(np.arange(Nslab, Natoms), i)
                        else:
                            index_other = np.delete(np.arange(Natoms), i)
                        blmin_array = np.array([self.blmin[num[i], num[j]] for j in index_other])
                        tf = np.min(atoms_mutated.get_distances(i, index_other,mic=True)/blmin_array) > 1.7
                    else:
                        tf = False
                    
                    if not tc and not tf:
                        break
                if tc or tf or tl:
                    atoms_mutated.positions[i] = posi_0
        return atoms_mutated


class RattleMutation_new(OffspringCreator):
    """ 
    A rattle mutation that that moves one atom and
    and places it in the vicinity of one of the other
    "top" atoms, randomly choosen.

    Parameters:

    blmin: Dictionary defining the minimum distance between atoms
    after the rattle.
    n_top: Number of atoms optimized by the GA.
    rattle_prop: The probability with which each atom is rattled.
    
    """
    def __init__(self, blmin, n_top, rattle_prop=0.4, min_z=None,
                 verbose=False, descriptor='RattleMutation_new'):
        OffspringCreator.__init__(self, verbose)
        self.blmin = blmin
        self.n_top = n_top
        self.rattle_prop = rattle_prop
        self.min_z = min_z
        self.descriptor = descriptor
        self.min_inputs = 1

    def get_new_individual(self, parents, **kwargs):
        f = parents[0]

        indi = self.mutate(f)
        if indi is None:
            return indi, 'mutation: rattle'
            
        indi = self.initialize_individual(f, indi)
        #indi.info['data']['parents'] = [f.info['confid']]  # remove for GP-search

        return self.finalize_individual(indi), 'mutation: rattle'

    def mutate(self, atoms):
        num = atoms.numbers
        atoms_mutated = atoms.copy()
        Natoms = len(atoms)
        Nslab = Natoms - self.n_top
        for i in range(Nslab, Natoms):
            if np.random.random() < self.rattle_prop:
                posi_0 = np.copy(atoms_mutated.positions[i])
                tc = False
                tf = False

                # Random order of atoms
                #atom_indicies = np.random.permutation(np.arange(Nslab, Natoms))
                atom_indicies = np.random.permutation(np.arange(Natoms))
                for j in atom_indicies:
                    if j==i:
                        continue
                    if self.min_z is not None:
                        if atoms_mutated.positions[j][2] < self.min_z:
                            continue
                    posi_j = np.copy(atoms_mutated.positions[j])
                    r_min = self.blmin[(num[i], num[j])]
                    r_max = 1.7*r_min
                    for k in range(20):
                        atoms_mutated.positions[i] = posi_j
                        # Then Rattle within a circle
                        
                        r = np.random.uniform(r_min**3, r_max**3)**(1/3)
                        theta = np.random.uniform(low=0, high=2*np.pi)
                        phi = np.random.uniform(low=0, high=np.pi)
                        pos_add = r * np.array([np.cos(theta)*np.sin(phi),
                                                np.sin(theta)*np.sin(phi),
                                                np.cos(phi)])
                        atoms_mutated.positions[i] += pos_add

                        # Too low - made to avoid that atom is placed below slab
                        if self.min_z is not None:
                            tl = atoms_mutated.positions[i][2] < self.min_z
                            if tl:
                                continue
                        else:
                            tl = False
                        
                        # too close
                        tc = atoms_too_close_two_sets(atoms_mutated[:i] + atoms_mutated[i+1:],
                                                      atoms_mutated[i], self.blmin)
                        # too far
                        index_other = np.delete(np.arange(Natoms), i)
                        blmin_array = np.array([self.blmin[num[i], num[j]] for j in index_other])
                        tf = np.min(atoms_mutated.get_distances(i, index_other, mic=True)/blmin_array) > 1.7
                        
                        if not tc and not tf:
                            break
                    if tc or tf or tl:
                        atoms_mutated.positions[i] = posi_0
        return atoms_mutated
    
class PermutationMutation(OffspringCreator):
    """Mutation that permutes a percentage of the atom types in the cluster.

       Parameters:

       n_top: Number of atoms optimized by the GA.
       probability: The probability with which an atom is permuted.
    """

    def __init__(self, blmin, n_top, probability=0.33, verbose=False, descriptor='PermutationMutation'):
        OffspringCreator.__init__(self, verbose)
        self.blmin = blmin
        self.n_top = n_top
        self.probability = probability
        self.descriptor = descriptor
        self.min_inputs = 1

    def get_new_individual(self, parents, **kwargs):
        f = parents[0]

        indi = self.mutate(f)
        if indi is None:
            return indi, 'mutation: permutation'
            
        indi = self.initialize_individual(f, indi)
        # indi.info['data']['parents'] = [f.info['confid']]  # remove for GP-search

        return self.finalize_individual(indi), 'mutation: permutation'

    def mutate(self, atoms):
        """ Does the actual mutation. """
        a = atoms.copy()
        Natoms = len(a)
        Nslab = Natoms - self.n_top
        pos = a.get_positions()[-self.n_top:]
        num = a.numbers
        num_top = num[-self.n_top:]
        num_unique = list(set(num_top))
        assert len(num_unique) > 1, 'Permutations with one atomic type is not valid'
        m = int(ceil(float(self.n_top) * self.probability / 2.))
        for _ in range(m):
            swap_succesfull = False
            for i in range(100):
                # Find index pair of atoms to swap
                i = j = 0
                while num[i] == num[j]:
                    i = np.random.randint(Nslab, Natoms)
                    j = np.random.randint(Nslab, Natoms)

                # Swap two atoms
                pos_i = a.positions[i].copy()
                pos_j = a.positions[j].copy()
                a.positions[i] = pos_j
                a.positions[j] = pos_i

                # Check if atoms are too close
                tc_i = atoms_too_close_two_sets(a[:i] + a[i+1:], a[i], self.blmin)
                tc_j = atoms_too_close_two_sets(a[:j] + a[j+1:], a[j], self.blmin)
                if tc_i or tc_j:
                    # reset swap
                    a.positions[i] = pos_i
                    a.positions[j] = pos_j
                    continue
                else:
                    swap_succesfull = True
                    break
            if not swap_succesfull:
                # swap_back
                a.positions[i] = pos_i
                a.positions[j] = pos_j
        return a


class MirrorMutation(OffspringCreator):
    """ A mirror mutation, as described in
        TO BE PUBLISHED.
        This mutation mirrors half of the cluster in a
        randomly oriented cutting plane discarding the other half.

        Parameters:
        blmin: Dictionary defining the minimum allowed
        distance between atoms.
        n_top: Number of atoms the GA optimizes.
        reflect: Defines if the mirrored half is also reflected
        perpendicular to the mirroring plane.

    """
    def __init__(self, blmin, n_top, reflect=False, verbose=False):
        OffspringCreator.__init__(self, verbose)
        self.blmin = blmin
        self.n_top = n_top
        self.reflect = reflect
        self.descriptor = 'MirrorMutation'
        self.min_inputs = 1

    def get_new_individual(self, parents):
        f = parents[0]

        indi = self.mutate(f)
        if indi is None:
            return indi, 'mutation: mirror'
            
        indi = self.initialize_individual(f, indi)
        indi.info['data']['parents'] = [f.info['confid']]

        return self.finalize_individual(indi), 'mutation: mirror'

    def mutate(self, atoms):
        """ Do the mutation of the atoms input. """

        reflect = self.reflect
        tc = True
        slab = atoms[0:len(atoms) - self.n_top]
        top = atoms[len(atoms) - self.n_top: len(atoms)]
        num = top.numbers
        unique_types = list(set(num))
        nu = dict()
        for u in unique_types:
            nu[u] = sum(num == u)
            
        n_tries = 1000
        counter = 0
        changed = False

        while tc and counter < n_tries:
            counter += 1
            cand = top.copy()
            pos = cand.get_positions()

            cm = np.average(top.get_positions(), axis=0)

            # first select a randomly oriented cutting plane
            theta = pi * np.random.random()
            phi = 2. * pi * np.random.random()
            n = (cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta))
            n = np.array(n)

            # Calculate all atoms signed distance to the cutting plane
            D = []
            for (i, p) in enumerate(pos):
                d = np.dot(p - cm, n)
                D.append((i, d))

            # Sort the atoms by their signed distance
            D.sort(key=lambda x: x[1])
            nu_taken = dict()

            # Select half of the atoms needed for a full cluster
            p_use = []
            n_use = []
            for (i, d) in D:
                if num[i] not in nu_taken.keys():
                    nu_taken[num[i]] = 0
                if nu_taken[num[i]] < nu[num[i]] / 2.:
                    p_use.append(pos[i])
                    n_use.append(num[i])
                    nu_taken[num[i]] += 1

            # calculate the mirrored position and add these.
            pn = []
            for p in p_use:
                pt = p - 2. * np.dot(p - cm, n) * n
                if reflect:
                    pt = -pt + 2 * cm + 2 * n * np.dot(pt - cm, n)
                pn.append(pt)

            n_use.extend(n_use)
            p_use.extend(pn)

            # In the case of an uneven number of
            # atoms we need to add one extra
            for n in nu.keys():
                if nu[n] % 2 == 0:
                    continue
                while sum(n_use == n) > nu[n]:
                    for i in range(int(len(n_use) / 2), len(n_use)):
                        if n_use[i] == n:
                            del p_use[i]
                            del n_use[i]
                            break
                assert sum(n_use == n) == nu[n]

            # Make sure we have the correct number of atoms
            # and rearrange the atoms so they are in the right order
            for i in range(len(n_use)):
                if num[i] == n_use[i]:
                    continue
                for j in range(i + 1, len(n_use)):
                    if n_use[j] == num[i]:
                        tn = n_use[i]
                        tp = p_use[i]
                        n_use[i] = n_use[j]
                        p_use[i] = p_use[j]
                        p_use[j] = tp
                        n_use[j] = tn

            # Finally we check that nothing is too close in the end product.
            cand = Atoms(num, p_use, cell=slab.get_cell(), pbc=slab.get_pbc())
            tc = atoms_too_close(cand, self.blmin)
            if tc:
                continue
            tc = atoms_too_close_two_sets(slab, cand, self.blmin)
            if not changed and counter > n_tries // 2:
                reflect = not reflect
                changed = True
            tot = slab + cand
        if counter == n_tries:
            return None
        return tot


if __name__ == '__main__':
    # Testing
    from ase.io import read
    from ase.visualize import view

    from ase.ga.utilities import closest_distances_generator
    from ase.ga.utilities import get_all_atom_types

    a = read('/home/mkb/DFT/gpLEA/gr_on_Ir/slab0/cand0.traj', index='0')

    unique_atom_numbers = list(set(a.get_atomic_numbers()))
    blmin = closest_distances_generator(atom_numbers = unique_atom_numbers,
                                       ratio_of_covalent_radii = 0.7)

    rattle = RattleMutation(blmin, n_top=35, rattle_strength=0.1, rattle_prop=0.1, rattle_mode='cylinder', cylinder_direction=0)

    a_rattle, _ = rattle.get_new_individual([a])
    view([a, a_rattle])