1from ase.ga.offspring_creator import OffspringCreator 

2from ase import Atoms 

3from itertools import chain 

4import numpy as np 

5 

6 

7class Crossover(OffspringCreator): 

8 """Base class for all particle crossovers. 

9 

10 Originally intended for medium sized particles 

11 

12 Do not call this class directly.""" 

13 

14 def __init__(self, rng=np.random): 

15 OffspringCreator.__init__(self, rng=rng) 

16 self.descriptor = 'Crossover' 

17 self.min_inputs = 2 

18 

19 

20class CutSpliceCrossover(Crossover): 

21 """Crossover that cuts two particles through a plane in space and 

22 merges two halfes from different particles together. 

23 

24 Implementation of the method presented in: 

25 D. M. Deaven and K. M. Ho, Phys. Rev. Lett., 75, 2, 288-291 (1995) 

26 

27 It keeps the correct composition by randomly assigning elements in 

28 the new particle. If some of the atoms in the two particle halves 

29 are too close, the halves are moved away from each other perpendicular 

30 to the cutting plane. 

31 

32 Parameters: 

33 

34 blmin: dictionary of minimum distance between atomic numbers. 

35 e.g. {(28,29): 1.5} 

36 

37 keep_composition: boolean that signifies if the composition should 

38 be the same as in the parents. 

39 

40 rng: Random number generator 

41 By default numpy.random. 

42 """ 

43 

44 def __init__(self, blmin, keep_composition=True, rng=np.random): 

45 Crossover.__init__(self, rng=rng) 

46 self.blmin = blmin 

47 self.keep_composition = keep_composition 

48 self.descriptor = 'CutSpliceCrossover' 

49 

50 def get_new_individual(self, parents): 

51 f, m = parents 

52 

53 indi = self.initialize_individual(f) 

54 indi.info['data']['parents'] = [i.info['confid'] for i in parents] 

55 

56 theta = self.rng.rand() * 2 * np.pi # 0,2pi 

57 phi = self.rng.rand() * np.pi # 0,pi 

58 e = np.array((np.sin(phi) * np.cos(theta), 

59 np.sin(theta) * np.sin(phi), 

60 np.cos(phi))) 

61 eps = 0.0001 

62 

63 f.translate(-f.get_center_of_mass()) 

64 m.translate(-m.get_center_of_mass()) 

65 

66 # Get the signed distance to the cutting plane 

67 # We want one side from f and the other side from m 

68 fmap = [np.dot(x, e) for x in f.get_positions()] 

69 mmap = [-np.dot(x, e) for x in m.get_positions()] 

70 ain = sorted([i for i in chain(fmap, mmap) if i > 0], 

71 reverse=True) 

72 aout = sorted([i for i in chain(fmap, mmap) if i < 0], 

73 reverse=True) 

74 

75 off = len(ain) - len(f) 

76 

77 # Translating f and m to get the correct number of atoms 

78 # in the offspring 

79 if off < 0: 

80 # too few 

81 # move f and m away from the plane 

82 dist = (abs(aout[abs(off) - 1]) + abs(aout[abs(off)])) * .5 

83 f.translate(e * dist) 

84 m.translate(-e * dist) 

85 elif off > 0: 

86 # too many 

87 # move f and m towards the plane 

88 dist = (abs(ain[-off - 1]) + abs(ain[-off])) * .5 

89 f.translate(-e * dist) 

90 m.translate(e * dist) 

91 if off != 0 and dist == 0: 

92 # Exactly same position => we continue with the wrong number 

93 # of atoms. What should be done? Fail or return None or 

94 # remove one of the two atoms with exactly the same position. 

95 pass 

96 

97 # Determine the contributing parts from f and m 

98 tmpf, tmpm = Atoms(), Atoms() 

99 for atom in f: 

100 if np.dot(atom.position, e) > 0: 

101 atom.tag = 1 

102 tmpf.append(atom) 

103 for atom in m: 

104 if np.dot(atom.position, e) < 0: 

105 atom.tag = 2 

106 tmpm.append(atom) 

107 

108 # Check that the correct composition is employed 

109 if self.keep_composition: 

110 opt_sm = sorted(f.numbers) 

111 tmpf_numbers = list(tmpf.numbers) 

112 tmpm_numbers = list(tmpm.numbers) 

113 cur_sm = sorted(tmpf_numbers + tmpm_numbers) 

114 # correct_by: dictionary that specifies how many 

115 # of the atom_numbers should be removed (a negative number) 

116 # or added (a positive number) 

117 correct_by = dict([(j, opt_sm.count(j)) for j in set(opt_sm)]) 

118 for n in cur_sm: 

119 correct_by[n] -= 1 

120 correct_in = tmpf if self.rng.choice([0, 1]) else tmpm 

121 to_add, to_rem = [], [] 

122 for num, amount in correct_by.items(): 

123 if amount > 0: 

124 to_add.extend([num] * amount) 

125 elif amount < 0: 

126 to_rem.extend([num] * abs(amount)) 

127 for add, rem in zip(to_add, to_rem): 

128 tbc = [a.index for a in correct_in if a.number == rem] 

129 if len(tbc) == 0: 

130 pass 

131 ai = self.rng.choice(tbc) 

132 correct_in[ai].number = add 

133 

134 # Move the contributing apart if any distance is below blmin 

135 maxl = 0. 

136 for sv, min_dist in self.get_vectors_below_min_dist(tmpf + tmpm): 

137 lsv = np.linalg.norm(sv) # length of shortest vector 

138 d = [-np.dot(e, sv)] * 2 

139 d[0] += np.sqrt(np.dot(e, sv)**2 - lsv**2 + min_dist**2) 

140 d[1] -= np.sqrt(np.dot(e, sv)**2 - lsv**2 + min_dist**2) 

141 l = sorted([abs(i) for i in d])[0] / 2. + eps 

142 if l > maxl: 

143 maxl = l 

144 tmpf.translate(e * maxl) 

145 tmpm.translate(-e * maxl) 

146 

147 # Put the two parts together 

148 for atom in chain(tmpf, tmpm): 

149 indi.append(atom) 

150 

151 parent_message = ':Parents {0} {1}'.format(f.info['confid'], 

152 m.info['confid']) 

153 return (self.finalize_individual(indi), 

154 self.descriptor + parent_message) 

155 

156 def get_numbers(self, atoms): 

157 """Returns the atomic numbers of the atoms object using only 

158 the elements defined in self.elements""" 

159 ac = atoms.copy() 

160 if self.elements is not None: 

161 del ac[[a.index for a in ac 

162 if a.symbol in self.elements]] 

163 return ac.numbers 

164 

165 def get_vectors_below_min_dist(self, atoms): 

166 """Generator function that returns each vector (between atoms) 

167 that is shorter than the minimum distance for those atom types 

168 (set during the initialization in blmin).""" 

169 norm = np.linalg.norm 

170 ap = atoms.get_positions() 

171 an = atoms.numbers 

172 for i in range(len(atoms)): 

173 pos = atoms[i].position 

174 for j, d in enumerate([norm(k - pos) for k in ap[i:]]): 

175 if d == 0: 

176 continue 

177 min_dist = self.blmin[tuple(sorted((an[i], an[j + i])))] 

178 if d < min_dist: 

179 yield atoms[i].position - atoms[j + i].position, min_dist 

180 

181 def get_shortest_dist_vector(self, atoms): 

182 norm = np.linalg.norm 

183 mind = 10000. 

184 ap = atoms.get_positions() 

185 for i in range(len(atoms)): 

186 pos = atoms[i].position 

187 for j, d in enumerate([norm(k - pos) for k in ap[i:]]): 

188 if d == 0: 

189 continue 

190 if d < mind: 

191 mind = d 

192 lowpair = (i, j + i) 

193 return atoms[lowpair[0]].position - atoms[lowpair[1]].position