1"""Soft-mutation operator and associated tools""" 

2import inspect 

3import json 

4import numpy as np 

5from ase.data import covalent_radii 

6from ase.neighborlist import NeighborList 

7from ase.ga.offspring_creator import OffspringCreator 

8from ase.ga.utilities import atoms_too_close, gather_atoms_by_tag 

9from scipy.spatial.distance import cdist 

10 

11 

12class TagFilter: 

13 """Filter which constrains same-tag atoms to behave 

14 like internally rigid moieties. 

15 """ 

16 def __init__(self, atoms): 

17 self.atoms = atoms 

18 gather_atoms_by_tag(self.atoms) 

19 self.tags = self.atoms.get_tags() 

20 self.unique_tags = np.unique(self.tags) 

21 self.n = len(self.unique_tags) 

22 

23 def get_positions(self): 

24 all_pos = self.atoms.get_positions() 

25 cop_pos = np.zeros((self.n, 3)) 

26 for i in range(self.n): 

27 indices = np.where(self.tags == self.unique_tags[i]) 

28 cop_pos[i] = np.average(all_pos[indices], axis=0) 

29 return cop_pos 

30 

31 def set_positions(self, positions, **kwargs): 

32 cop_pos = self.get_positions() 

33 all_pos = self.atoms.get_positions() 

34 assert np.all(np.shape(positions) == np.shape(cop_pos)) 

35 for i in range(self.n): 

36 indices = np.where(self.tags == self.unique_tags[i]) 

37 shift = positions[i] - cop_pos[i] 

38 all_pos[indices] += shift 

39 self.atoms.set_positions(all_pos, **kwargs) 

40 

41 def get_forces(self, *args, **kwargs): 

42 f = self.atoms.get_forces() 

43 forces = np.zeros((self.n, 3)) 

44 for i in range(self.n): 

45 indices = np.where(self.tags == self.unique_tags[i]) 

46 forces[i] = np.sum(f[indices], axis=0) 

47 return forces 

48 

49 def get_masses(self): 

50 m = self.atoms.get_masses() 

51 masses = np.zeros(self.n) 

52 for i in range(self.n): 

53 indices = np.where(self.tags == self.unique_tags[i]) 

54 masses[i] = np.sum(m[indices]) 

55 return masses 

56 

57 def __len__(self): 

58 return self.n 

59 

60 

61class PairwiseHarmonicPotential: 

62 """Parent class for interatomic potentials of the type 

63 E(r_ij) = 0.5 * k_ij * (r_ij - r0_ij) ** 2 

64 """ 

65 def __init__(self, atoms, rcut=10.): 

66 self.atoms = atoms 

67 self.pos0 = atoms.get_positions() 

68 self.rcut = rcut 

69 

70 # build neighborlist 

71 nat = len(self.atoms) 

72 self.nl = NeighborList([self.rcut / 2.] * nat, skin=0., bothways=True, 

73 self_interaction=False) 

74 self.nl.update(self.atoms) 

75 

76 self.calculate_force_constants() 

77 

78 def calculate_force_constants(self): 

79 msg = 'Child class needs to define a calculate_force_constants() ' \ 

80 'method which computes the force constants and stores them ' \ 

81 'in self.force_constants (as a list which contains, for every ' \ 

82 'atom, a list of the atom\'s force constants with its neighbors.' 

83 raise NotImplementedError(msg) 

84 

85 def get_forces(self, atoms): 

86 pos = atoms.get_positions() 

87 cell = atoms.get_cell() 

88 forces = np.zeros_like(pos) 

89 

90 for i, p in enumerate(pos): 

91 indices, offsets = self.nl.get_neighbors(i) 

92 p = pos[indices] + np.dot(offsets, cell) 

93 r = cdist(p, [pos[i]]) 

94 v = (p - pos[i]) / r 

95 p0 = self.pos0[indices] + np.dot(offsets, cell) 

96 r0 = cdist(p0, [self.pos0[i]]) 

97 dr = r - r0 

98 forces[i] = np.dot(self.force_constants[i].T, dr * v) 

99 

100 return forces 

101 

102 

103def get_number_of_valence_electrons(Z): 

104 """Return the number of valence electrons for the element with 

105 atomic number Z, simply based on its periodic table group. 

106 """ 

107 groups = [[], [1, 3, 11, 19, 37, 55, 87], [2, 4, 12, 20, 38, 56, 88], 

108 [21, 39, 57, 89]] 

109 

110 for i in range(9): 

111 groups.append(i + np.array([22, 40, 72, 104])) 

112 

113 for i in range(6): 

114 groups.append(i + np.array([5, 13, 31, 49, 81, 113])) 

115 

116 for i, group in enumerate(groups): 

117 if Z in group: 

118 nval = i if i < 13 else i - 10 

119 break 

120 else: 

121 raise ValueError('Z=%d not included in this dataset.' % Z) 

122 

123 return nval 

124 

125 

126class BondElectroNegativityModel(PairwiseHarmonicPotential): 

127 """Pairwise harmonic potential where the force constants are 

128 determined using the "bond electronegativity" model, see: 

129 

130 * `Lyakhov, Oganov, Valle, Comp. Phys. Comm. 181 (2010) 1623-1632`__ 

131 

132 __ https://dx.doi.org/10.1016/j.cpc.2010.06.007 

133 

134 * `Lyakhov, Oganov, Phys. Rev. B 84 (2011) 092103`__ 

135 

136 __ https://dx.doi.org/10.1103/PhysRevB.84.092103 

137 """ 

138 def calculate_force_constants(self): 

139 cell = self.atoms.get_cell() 

140 pos = self.atoms.get_positions() 

141 num = self.atoms.get_atomic_numbers() 

142 nat = len(self.atoms) 

143 nl = self.nl 

144 

145 # computing the force constants 

146 s_norms = [] 

147 valence_states = [] 

148 r_cov = [] 

149 for i in range(nat): 

150 indices, offsets = nl.get_neighbors(i) 

151 p = pos[indices] + np.dot(offsets, cell) 

152 r = cdist(p, [pos[i]]) 

153 r_ci = covalent_radii[num[i]] 

154 s = 0. 

155 for j, index in enumerate(indices): 

156 d = r[j] - r_ci - covalent_radii[num[index]] 

157 s += np.exp(-d / 0.37) 

158 s_norms.append(s) 

159 valence_states.append(get_number_of_valence_electrons(num[i])) 

160 r_cov.append(r_ci) 

161 

162 self.force_constants = [] 

163 for i in range(nat): 

164 indices, offsets = nl.get_neighbors(i) 

165 p = pos[indices] + np.dot(offsets, cell) 

166 r = cdist(p, [pos[i]])[:, 0] 

167 fc = [] 

168 for j, ii in enumerate(indices): 

169 d = r[j] - r_cov[i] - r_cov[ii] 

170 chi_ik = 0.481 * valence_states[i] / (r_cov[i] + 0.5 * d) 

171 chi_jk = 0.481 * valence_states[ii] / (r_cov[ii] + 0.5 * d) 

172 cn_ik = s_norms[i] / np.exp(-d / 0.37) 

173 cn_jk = s_norms[ii] / np.exp(-d / 0.37) 

174 fc.append(np.sqrt(chi_ik * chi_jk / (cn_ik * cn_jk))) 

175 self.force_constants.append(np.array(fc)) 

176 

177 

178class SoftMutation(OffspringCreator): 

179 """Mutates the structure by displacing it along the lowest 

180 (nonzero) frequency modes found by vibrational analysis, as in: 

181 

182 `Lyakhov, Oganov, Valle, Comp. Phys. Comm. 181 (2010) 1623-1632`__ 

183 

184 __ https://dx.doi.org/10.1016/j.cpc.2010.06.007 

185 

186 As in the reference above, the next-lowest mode is used if the 

187 structure has already been softmutated along the current-lowest 

188 mode. This mutation hence acts in a deterministic way, in contrast 

189 to most other genetic operators. 

190 

191 If you find this implementation useful in your work, 

192 please consider citing: 

193 

194 `Van den Bossche, Gronbeck, Hammer, J. Chem. Theory Comput. 14 (2018)`__ 

195 

196 __ https://dx.doi.org/10.1021/acs.jctc.8b00039 

197 

198 in addition to the paper mentioned above. 

199 

200 Parameters: 

201 

202 blmin: dict 

203 The closest allowed interatomic distances on the form: 

204 {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers. 

205 

206 bounds: list 

207 Lower and upper limits (in Angstrom) for the largest 

208 atomic displacement in the structure. For a given mode, 

209 the algorithm starts at zero amplitude and increases 

210 it until either blmin is violated or the largest 

211 displacement exceeds the provided upper bound). 

212 If the largest displacement in the resulting structure 

213 is lower than the provided lower bound, the mutant is 

214 considered too similar to the parent and None is 

215 returned. 

216 

217 calculator: ASE calculator object 

218 The calculator to be used in the vibrational 

219 analysis. The default is to use a calculator 

220 based on pairwise harmonic potentials with force 

221 constants from the "bond electronegativity" 

222 model described in the reference above. 

223 Any calculator with a working :func:`get_forces()` 

224 method will work. 

225 

226 rcut: float 

227 Cutoff radius in Angstrom for the pairwise harmonic 

228 potentials. 

229 

230 used_modes_file: str or None 

231 Name of json dump file where previously used 

232 modes will be stored (and read). If None, 

233 no such file will be used. Default is to use 

234 the filename 'used_modes.json'. 

235 

236 use_tags: boolean 

237 Whether to use the atomic tags to preserve molecular identity. 

238 """ 

239 def __init__(self, blmin, bounds=[0.5, 2.0], 

240 calculator=BondElectroNegativityModel, rcut=10., 

241 used_modes_file='used_modes.json', use_tags=False, 

242 verbose=False): 

243 OffspringCreator.__init__(self, verbose) 

244 self.blmin = blmin 

245 self.bounds = bounds 

246 self.calc = calculator 

247 self.rcut = rcut 

248 self.used_modes_file = used_modes_file 

249 self.use_tags = use_tags 

250 self.descriptor = 'SoftMutation' 

251 

252 self.used_modes = {} 

253 if self.used_modes_file is not None: 

254 try: 

255 self.read_used_modes(self.used_modes_file) 

256 except IOError: 

257 # file doesn't exist (yet) 

258 pass 

259 

260 def _get_hessian(self, atoms, dx): 

261 """Returns the Hessian matrix d2E/dxi/dxj using a first-order 

262 central difference scheme with displacements dx. 

263 """ 

264 N = len(atoms) 

265 pos = atoms.get_positions() 

266 hessian = np.zeros((3 * N, 3 * N)) 

267 

268 for i in range(3 * N): 

269 row = np.zeros(3 * N) 

270 for direction in [-1, 1]: 

271 disp = np.zeros(3) 

272 disp[i % 3] = direction * dx 

273 pos_disp = np.copy(pos) 

274 pos_disp[i // 3] += disp 

275 atoms.set_positions(pos_disp) 

276 f = atoms.get_forces() 

277 row += -1 * direction * f.flatten() 

278 

279 row /= (2. * dx) 

280 hessian[i] = row 

281 

282 hessian += np.copy(hessian).T 

283 hessian *= 0.5 

284 atoms.set_positions(pos) 

285 

286 return hessian 

287 

288 def _calculate_normal_modes(self, atoms, dx=0.02, massweighing=False): 

289 """Performs the vibrational analysis.""" 

290 hessian = self._get_hessian(atoms, dx) 

291 if massweighing: 

292 m = np.array([np.repeat(atoms.get_masses()**-0.5, 3)]) 

293 hessian *= (m * m.T) 

294 

295 eigvals, eigvecs = np.linalg.eigh(hessian) 

296 modes = {eigval: eigvecs[:, i] for i, eigval in enumerate(eigvals)} 

297 return modes 

298 

299 def animate_mode(self, atoms, mode, nim=30, amplitude=1.0): 

300 """Returns an Atoms object showing an animation of the mode.""" 

301 pos = atoms.get_positions() 

302 mode = mode.reshape(np.shape(pos)) 

303 animation = [] 

304 for i in range(nim): 

305 newpos = pos + amplitude * mode * np.sin(i * 2 * np.pi / nim) 

306 image = atoms.copy() 

307 image.positions = newpos 

308 animation.append(image) 

309 return animation 

310 

311 def read_used_modes(self, filename): 

312 """Read used modes from json file.""" 

313 with open(filename, 'r') as fd: 

314 modes = json.load(fd) 

315 self.used_modes = {int(k): modes[k] for k in modes} 

316 return 

317 

318 def write_used_modes(self, filename): 

319 """Dump used modes to json file.""" 

320 with open(filename, 'w') as fd: 

321 json.dump(self.used_modes, fd) 

322 return 

323 

324 def get_new_individual(self, parents): 

325 f = parents[0] 

326 

327 indi = self.mutate(f) 

328 if indi is None: 

329 return indi, 'mutation: soft' 

330 

331 indi = self.initialize_individual(f, indi) 

332 indi.info['data']['parents'] = [f.info['confid']] 

333 

334 return self.finalize_individual(indi), 'mutation: soft' 

335 

336 def mutate(self, atoms): 

337 """Does the actual mutation.""" 

338 a = atoms.copy() 

339 

340 if inspect.isclass(self.calc): 

341 assert issubclass(self.calc, PairwiseHarmonicPotential) 

342 calc = self.calc(atoms, rcut=self.rcut) 

343 else: 

344 calc = self.calc 

345 a.calc = calc 

346 

347 if self.use_tags: 

348 a = TagFilter(a) 

349 

350 pos = a.get_positions() 

351 modes = self._calculate_normal_modes(a) 

352 

353 # Select the mode along which we want to move the atoms; 

354 # The first 3 translational modes as well as previously 

355 # applied modes are discarded. 

356 

357 keys = np.array(sorted(modes)) 

358 index = 3 

359 confid = atoms.info['confid'] 

360 if confid in self.used_modes: 

361 while index in self.used_modes[confid]: 

362 index += 1 

363 self.used_modes[confid].append(index) 

364 else: 

365 self.used_modes[confid] = [index] 

366 

367 if self.used_modes_file is not None: 

368 self.write_used_modes(self.used_modes_file) 

369 

370 key = keys[index] 

371 mode = modes[key].reshape(np.shape(pos)) 

372 

373 # Find a suitable amplitude for translation along the mode; 

374 # at every trial amplitude both positive and negative 

375 # directions are tried. 

376 

377 mutant = atoms.copy() 

378 amplitude = 0. 

379 increment = 0.1 

380 direction = 1 

381 largest_norm = np.max(np.apply_along_axis(np.linalg.norm, 1, mode)) 

382 

383 def expand(atoms, positions): 

384 if isinstance(atoms, TagFilter): 

385 a.set_positions(positions) 

386 return a.atoms.get_positions() 

387 else: 

388 return positions 

389 

390 while amplitude * largest_norm < self.bounds[1]: 

391 pos_new = pos + direction * amplitude * mode 

392 pos_new = expand(a, pos_new) 

393 mutant.set_positions(pos_new) 

394 mutant.wrap() 

395 too_close = atoms_too_close(mutant, self.blmin, 

396 use_tags=self.use_tags) 

397 if too_close: 

398 amplitude -= increment 

399 pos_new = pos + direction * amplitude * mode 

400 pos_new = expand(a, pos_new) 

401 mutant.set_positions(pos_new) 

402 mutant.wrap() 

403 break 

404 

405 if direction == 1: 

406 direction = -1 

407 else: 

408 direction = 1 

409 amplitude += increment 

410 

411 if amplitude * largest_norm < self.bounds[0]: 

412 mutant = None 

413 

414 return mutant