candidate_generation.py~ 3.98 KB
Newer Older
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
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
import numpy as np
from abc import ABC, abstractmethod
from ase.data import covalent_radii
from ase.geometry import get_distances

class OffspringCreator(ABC):
    def __init__(self, blmin=0.7, blmax=1.3):
        self.blmin = blmin
        self.blmax = blmax
        self.description = 'Unspecified'

    @abstractmethod
    def get_new_candidate(self):
        pass

    def check_valid_bondlengths(self, a, indices=None):
        bl = self.get_distances_as_fraction_of_covalent(a,indices)

        # Filter away self interactions.
        bl = bl[bl>1e-6].reshape(bl.shape[0],bl.shape[1]-1)

        # Check if atoms are too close
        tc = np.any(bl < self.blmin)

        # Check if there are isolated atoms
        if self.blmax is not None:
            isolated = np.any(np.all(bl > self.blmax, axis=1))
        else:
            isolated = False

        is_valid = not tc and not isolated
        return is_valid

    def get_covalent_distance_from_atom_numbers(self, a, indices=None):
        r_cov = np.array([covalent_radii[n] for n in a.get_atomic_numbers()])
        if indices is None:
            r_cov_sub = r_cov
        else:
            r_cov_sub = r_cov[indices]
        cd_mat = r_cov_sub.reshape(-1,1) + r_cov.reshape(1,-1)
        return cd_mat
    
    def get_distances_as_fraction_of_covalent(self, a, indices=None, covalent_distances=None):
        if indices is None:
            indices = np.arange(len(a))

        if covalent_distances is None:
            cd = self.get_covalent_distance_from_atom_numbers(a, indices=indices)
        else:
            cd = covalent_distances[indices,:]
        _, d = get_distances(a[indices].positions,
                             a.positions,
                             cell=a.get_cell(),
                             pbc=a.get_pbc())
        bl = d/cd
        
        return bl

    def finalize(self, a):
        a.info['key_value_pairs'] = {'origin': self.description}
        valid_bondlengths = self.check_valid_bondlengths(a)
        assert valid_bondlengths
        return a

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
    
class RattleMutation(OffspringCreator):
    def __init__(self, n_top, Nrattle=3, rattle_range=3, blmin=0.7, blmax=1.3):
        OffspringCreator.__init__(self, blmin=blmin, blmax=blmax)
        self.description = 'RattleMutation'
        self.n_top = n_top
        self.probability = Nrattle/n_top
        self.rattle_range = rattle_range

    def get_new_candidate(self, parents):
        a = parents[0]
        a = self.rattle(a)
        a = self.finalize(a)
        return a

    def rattle(self, atoms):
        a = atoms.copy()
        Natoms = len(a)
        Nslab = Natoms - self.n_top

        for i in range(Nslab,Natoms):
            if np.random.random() < self.probability:
                for _ in range(100): 
                    posi_0 = np.copy(a.positions[i])

                    # Perform rattle
                    pos_add = pos_add_sphere(self.rattle_range)
                    a.positions[i] += pos_add
                    
                    # Check if rattle was valid
                    valid_bondlengths = self.check_valid_bondlengths(a, indices=[i])

                    if not valid_bondlengths:
                        a.positions[i] = posi_0
                    else:
                        break
        return a
                    

class candidateOperation():
    def __init__(self):
        pass


if __name__ == '__main__':
    from ase.io import read
    from ase.visualize import view
    
    a = read('/home/mkb/DFT/gpLEA/Si3x3/ref/gm_unrelaxed_done.traj', index='0')
    
    rattle = RattleMutation(n_top=16, Nrattle=3)

    a_mut = rattle.get_new_candidate([a])
    view([a,a_mut])