Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""Various utility methods used troughout the GA.""" 

2import os 

3import time 

4import itertools 

5import numpy as np 

6from scipy.spatial.distance import cdist 

7from ase.io import write, read 

8from ase.geometry.cell import cell_to_cellpar 

9from ase.geometry.rdf import get_rdf 

10from ase.data import covalent_radii 

11from ase.ga import get_neighbor_list 

12 

13 

14def closest_distances_generator(atom_numbers, ratio_of_covalent_radii): 

15 """Generates the blmin dict used across the GA. 

16 The distances are based on the covalent radii of the atoms. 

17 """ 

18 cr = covalent_radii 

19 ratio = ratio_of_covalent_radii 

20 

21 blmin = dict() 

22 for i in atom_numbers: 

23 blmin[(i, i)] = cr[i] * 2 * ratio 

24 for j in atom_numbers: 

25 if i == j: 

26 continue 

27 if (i, j) in blmin.keys(): 

28 continue 

29 blmin[(i, j)] = blmin[(j, i)] = ratio * (cr[i] + cr[j]) 

30 return blmin 

31 

32 

33def get_mic_distance(p1, p2, cell, pbc): 

34 """This method calculates the shortest distance between p1 and p2 

35 through the cell boundaries defined by cell and pbc. 

36 This method works for reasonable unit cells, but not for extremely 

37 elongated ones. 

38 """ 

39 ct = cell.T 

40 pos = np.array((p1, p2)) 

41 scaled = np.linalg.solve(ct, pos.T).T 

42 for i in range(3): 

43 if pbc[i]: 

44 scaled[:, i] %= 1.0 

45 scaled[:, i] %= 1.0 

46 P = np.dot(scaled, cell) 

47 

48 pbc_directions = [[-1, 1] * int(direction) + [0] for direction in pbc] 

49 translations = np.array(list(itertools.product(*pbc_directions))).T 

50 p0r = np.tile(np.reshape(P[0, :], (3, 1)), (1, translations.shape[1])) 

51 p1r = np.tile(np.reshape(P[1, :], (3, 1)), (1, translations.shape[1])) 

52 dp_vec = p0r + np.dot(ct, translations) 

53 d = np.min(np.power(p1r - dp_vec, 2).sum(axis=0))**0.5 

54 return d 

55 

56 

57def db_call_with_error_tol(db_cursor, expression, args=[]): 

58 """In case the GA is used on older versions of networking 

59 filesystems there might be some delays. For this reason 

60 some extra error tolerance when calling the SQLite db is 

61 employed. 

62 """ 

63 import sqlite3 

64 i = 0 

65 while i < 10: 

66 try: 

67 db_cursor.execute(expression, args) 

68 return 

69 except sqlite3.OperationalError as e: 

70 print(e) 

71 time.sleep(2.) 

72 i += 1 

73 raise sqlite3.OperationalError( 

74 'Database still locked after 10 attempts (20 s)') 

75 

76 

77def save_trajectory(confid, trajectory, folder): 

78 """Saves traj files to the database folder. 

79 This method should never be used directly, 

80 but only through the DataConnection object. 

81 """ 

82 fname = os.path.join(folder, 'traj%05d.traj' % confid) 

83 write(fname, trajectory) 

84 return fname 

85 

86 

87def get_trajectory(fname): 

88 """Extra error tolerance when loading traj files.""" 

89 fname = str(fname) 

90 try: 

91 t = read(fname) 

92 except IOError as e: 

93 print('get_trajectory error ' + e) 

94 return t 

95 

96 

97def gather_atoms_by_tag(atoms): 

98 """Translates same-tag atoms so that they lie 'together', 

99 with distance vectors as in the minimum image convention.""" 

100 tags = atoms.get_tags() 

101 pos = atoms.get_positions() 

102 for tag in list(set(tags)): 

103 indices = np.where(tags == tag)[0] 

104 if len(indices) == 1: 

105 continue 

106 vectors = atoms.get_distances(indices[0], indices[1:], 

107 mic=True, vector=True) 

108 pos[indices[1:]] = pos[indices[0]] + vectors 

109 atoms.set_positions(pos) 

110 

111 

112def atoms_too_close(atoms, bl, use_tags=False): 

113 """Checks if any atoms in a are too close, as defined by 

114 the distances in the bl dictionary. 

115 

116 use_tags: whether to use the Atoms tags to disable distance 

117 checking within a set of atoms with the same tag. 

118 

119 Note: if certain atoms are constrained and use_tags is True, 

120 this method may return unexpected results in case the 

121 contraints prevent same-tag atoms to be gathered together in 

122 the minimum-image-convention. In such cases, one should 

123 (1) release the relevant constraints, 

124 (2) apply the gather_atoms_by_tag function, and 

125 (3) re-apply the constraints, before using the 

126 atoms_too_close function. 

127 """ 

128 a = atoms.copy() 

129 if use_tags: 

130 gather_atoms_by_tag(a) 

131 

132 pbc = a.get_pbc() 

133 cell = a.get_cell() 

134 num = a.get_atomic_numbers() 

135 pos = a.get_positions() 

136 tags = a.get_tags() 

137 unique_types = sorted(list(set(num))) 

138 

139 neighbours = [] 

140 for i in range(3): 

141 if pbc[i]: 

142 neighbours.append([-1, 0, 1]) 

143 else: 

144 neighbours.append([0]) 

145 

146 for nx, ny, nz in itertools.product(*neighbours): 

147 displacement = np.dot(cell.T, np.array([nx, ny, nz]).T) 

148 pos_new = pos + displacement 

149 distances = cdist(pos, pos_new) 

150 

151 if nx == 0 and ny == 0 and nz == 0: 

152 if use_tags and len(a) > 1: 

153 x = np.array([tags]).T 

154 distances += 1e2 * (cdist(x, x) == 0) 

155 else: 

156 distances += 1e2 * np.identity(len(a)) 

157 

158 iterator = itertools.combinations_with_replacement(unique_types, 2) 

159 for type1, type2 in iterator: 

160 x1 = np.where(num == type1) 

161 x2 = np.where(num == type2) 

162 if np.min(distances[x1].T[x2]) < bl[(type1, type2)]: 

163 return True 

164 

165 return False 

166 

167 

168def atoms_too_close_two_sets(a, b, bl): 

169 """Checks if any atoms in a are too close to an atom in b, 

170 as defined by the bl dictionary.""" 

171 pbc_a = a.get_pbc() 

172 pbc_b = b.get_pbc() 

173 cell_a = a.get_cell() 

174 cell_b = a.get_cell() 

175 assert np.allclose(pbc_a, pbc_b), (pbc_a, pbc_b) 

176 assert np.allclose(cell_a, cell_b), (cell_a, cell_b) 

177 

178 pos_a = a.get_positions() 

179 pos_b = b.get_positions() 

180 

181 num_a = a.get_atomic_numbers() 

182 num_b = b.get_atomic_numbers() 

183 unique_types = sorted(set(list(num_a) + list(num_b))) 

184 

185 neighbours = [] 

186 for i in range(3): 

187 neighbours.append([-1, 0, 1] if pbc_a[i] else [0]) 

188 

189 for nx, ny, nz in itertools.product(*neighbours): 

190 displacement = np.dot(cell_a.T, np.array([nx, ny, nz]).T) 

191 pos_b_disp = pos_b + displacement 

192 distances = cdist(pos_a, pos_b_disp) 

193 

194 for type1 in unique_types: 

195 if type1 not in num_a: 

196 continue 

197 x1 = np.where(num_a == type1) 

198 

199 for type2 in unique_types: 

200 if type2 not in num_b: 

201 continue 

202 x2 = np.where(num_b == type2) 

203 if np.min(distances[x1].T[x2]) < bl[(type1, type2)]: 

204 return True 

205 return False 

206 

207 

208def get_all_atom_types(slab, atom_numbers_to_optimize): 

209 """Utility method used to extract all unique atom types 

210 from the atoms object slab and the list of atomic numbers 

211 atom_numbers_to_optimize. 

212 """ 

213 from_slab = list(set(slab.numbers)) 

214 from_top = list(set(atom_numbers_to_optimize)) 

215 from_slab.extend(from_top) 

216 return list(set(from_slab)) 

217 

218 

219def get_distance_matrix(atoms, self_distance=1000): 

220 """NB: This function is way slower than atoms.get_all_distances() 

221 

222 Returns a numpy matrix with the distances between the atoms 

223 in the supplied atoms object, with the indices of the matrix 

224 corresponding to the indices in the atoms object. 

225 

226 The parameter self_distance will be put in the diagonal 

227 elements ([i][i]) 

228 """ 

229 dm = np.zeros([len(atoms), len(atoms)]) 

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

231 dm[i][i] = self_distance 

232 for j in range(i + 1, len(atoms)): 

233 rij = atoms.get_distance(i, j) 

234 dm[i][j] = rij 

235 dm[j][i] = rij 

236 return dm 

237 

238 

239def get_nndist(atoms, distance_matrix): 

240 """Returns an estimate of the nearest neighbor bond distance 

241 in the supplied atoms object given the supplied distance_matrix. 

242 

243 The estimate comes from the first peak in the radial distribution 

244 function. 

245 """ 

246 rmax = 10. # No bonds longer than 10 angstrom expected 

247 nbins = 200 

248 rdf, dists = get_rdf(atoms, rmax, nbins, distance_matrix) 

249 return dists[np.argmax(rdf)] 

250 

251 

252def get_nnmat(atoms, mic=False): 

253 """Calculate the nearest neighbor matrix as specified in 

254 S. Lysgaard et al., Top. Catal., 2014, 57 (1-4), pp 33-39 

255 

256 Returns an array of average numbers of nearest neighbors 

257 the order is determined by self.elements. 

258 Example: self.elements = ["Cu", "Ni"] 

259 get_nnmat returns a single list [Cu-Cu bonds/N(Cu), 

260 Cu-Ni bonds/N(Cu), Ni-Cu bonds/N(Ni), Ni-Ni bonds/N(Ni)] 

261 where N(element) is the number of atoms of the type element 

262 in the atoms object. 

263 

264 The distance matrix can be quite costly to calculate every 

265 time nnmat is required (and disk intensive if saved), thus 

266 it makes sense to calculate nnmat along with e.g. the 

267 potential energy and save it in atoms.info['data']['nnmat']. 

268 """ 

269 if 'data' in atoms.info and 'nnmat' in atoms.info['data']: 

270 return atoms.info['data']['nnmat'] 

271 elements = sorted(set(atoms.get_chemical_symbols())) 

272 nnmat = np.zeros((len(elements), len(elements))) 

273 # dm = get_distance_matrix(atoms) 

274 dm = atoms.get_all_distances(mic=mic) 

275 nndist = get_nndist(atoms, dm) + 0.2 

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

277 row = [j for j in range(len(elements)) 

278 if atoms[i].symbol == elements[j]][0] 

279 neighbors = [j for j in range(len(dm[i])) if dm[i][j] < nndist] 

280 for n in neighbors: 

281 column = [j for j in range(len(elements)) 

282 if atoms[n].symbol == elements[j]][0] 

283 nnmat[row][column] += 1 

284 # divide by the number of that type of atoms in the structure 

285 for i, el in enumerate(elements): 

286 nnmat[i] /= len([j for j in range(len(atoms)) 

287 if atoms[int(j)].symbol == el]) 

288 # makes a single list out of a list of lists 

289 nnlist = np.reshape(nnmat, (len(nnmat)**2)) 

290 return nnlist 

291 

292 

293def get_nnmat_string(atoms, decimals=2, mic=False): 

294 nnmat = get_nnmat(atoms, mic=mic) 

295 s = '-'.join(['{1:2.{0}f}'.format(decimals, i) 

296 for i in nnmat]) 

297 if len(nnmat) == 1: 

298 return s + '-' 

299 return s 

300 

301 

302def get_connections_index(atoms, max_conn=5, no_count_types=None): 

303 """This method returns a dictionary where each key value are a 

304 specific number of neighbors and list of atoms indices with 

305 that amount of neighbors respectively. The method utilizes the 

306 neighbor list and hence inherit the restrictions for 

307 neighbors. Option added to remove connections between 

308 defined atom types. 

309 

310 Parameters 

311 ---------- 

312 

313 atoms : Atoms object 

314 The connections will be counted using this supplied Atoms object 

315 

316 max_conn : int 

317 Any atom with more connections than this will be counted as 

318 having max_conn connections. 

319 Default 5 

320 

321 no_count_types : list or None 

322 List of atomic numbers that should be excluded in the count. 

323 Default None (meaning all atoms count). 

324 """ 

325 conn = get_neighbor_list(atoms) 

326 

327 if conn is None: 

328 conn = get_neighborlist(atoms) 

329 

330 if no_count_types is None: 

331 no_count_types = [] 

332 

333 conn_index = {} 

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

335 if atoms[i].number not in no_count_types: 

336 cconn = min(len(conn[i]), max_conn - 1) 

337 if cconn not in conn_index: 

338 conn_index[cconn] = [] 

339 conn_index[cconn].append(i) 

340 

341 return conn_index 

342 

343 

344def get_atoms_connections(atoms, max_conn=5, no_count_types=None): 

345 """This method returns a list of the numbers of atoms 

346 with X number of neighbors. The method utilizes the 

347 neighbor list and hence inherit the restrictions for 

348 neighbors. Option added to remove connections between 

349 defined atom types. 

350 """ 

351 conn_index = get_connections_index(atoms, max_conn=max_conn, 

352 no_count_types=no_count_types) 

353 

354 no_of_conn = [0] * max_conn 

355 for i in conn_index: 

356 no_of_conn[i] += len(conn_index[i]) 

357 

358 return no_of_conn 

359 

360 

361def get_angles_distribution(atoms, ang_grid=9): 

362 """Method to get the distribution of bond angles 

363 in bins (default 9) with bonds defined from 

364 the get_neighbor_list(). 

365 """ 

366 conn = get_neighbor_list(atoms) 

367 

368 if conn is None: 

369 conn = get_neighborlist(atoms) 

370 

371 bins = [0] * ang_grid 

372 

373 for atom in atoms: 

374 for i in conn[atom.index]: 

375 for j in conn[atom.index]: 

376 if j != i: 

377 a = atoms.get_angle(i, atom.index, j) 

378 for k in range(ang_grid): 

379 if (k + 1) * 180. / ang_grid > a > k * 180. / ang_grid: 

380 bins[k] += 1 

381 # Removing dobbelt counting 

382 for i in range(ang_grid): 

383 bins[i] /= 2. 

384 return bins 

385 

386 

387def get_neighborlist(atoms, dx=0.2, no_count_types=None): 

388 """Method to get the a dict with list of neighboring 

389 atoms defined as the two covalent radii + fixed distance. 

390 Option added to remove neighbors between defined atom types. 

391 """ 

392 cell = atoms.get_cell() 

393 pbc = atoms.get_pbc() 

394 

395 if no_count_types is None: 

396 no_count_types = [] 

397 

398 conn = {} 

399 for atomi in atoms: 

400 conn_this_atom = [] 

401 for atomj in atoms: 

402 if atomi.index != atomj.index: 

403 if atomi.number not in no_count_types: 

404 if atomj.number not in no_count_types: 

405 d = get_mic_distance(atomi.position, 

406 atomj.position, 

407 cell, 

408 pbc) 

409 cri = covalent_radii[atomi.number] 

410 crj = covalent_radii[atomj.number] 

411 d_max = crj + cri + dx 

412 if d < d_max: 

413 conn_this_atom.append(atomj.index) 

414 conn[atomi.index] = conn_this_atom 

415 return conn 

416 

417 

418def get_atoms_distribution(atoms, number_of_bins=5, max_distance=8, 

419 center=None, no_count_types=None): 

420 """Method to get the distribution of atoms in the 

421 structure in bins of distances from a defined 

422 center. Option added to remove counting of 

423 certain atom types. 

424 """ 

425 pbc = atoms.get_pbc() 

426 cell = atoms.get_cell() 

427 if center is None: 

428 # Center used for the atom distribution if None is supplied! 

429 cx = sum(cell[:, 0]) / 2. 

430 cy = sum(cell[:, 1]) / 2. 

431 cz = sum(cell[:, 2]) / 2. 

432 center = (cx, cy, cz) 

433 bins = [0] * number_of_bins 

434 

435 if no_count_types is None: 

436 no_count_types = [] 

437 

438 for atom in atoms: 

439 if atom.number not in no_count_types: 

440 d = get_mic_distance(atom.position, center, cell, pbc) 

441 for k in range(number_of_bins - 1): 

442 min_dis_cur_bin = k * max_distance / (number_of_bins - 1.) 

443 max_dis_cur_bin = ((k + 1) * max_distance / 

444 (number_of_bins - 1.)) 

445 if min_dis_cur_bin < d < max_dis_cur_bin: 

446 bins[k] += 1 

447 if d > max_distance: 

448 bins[number_of_bins - 1] += 1 

449 return bins 

450 

451 

452def get_rings(atoms, rings=[5, 6, 7]): 

453 """This method return a list of the number of atoms involved 

454 in rings in the structures. It uses the neighbor 

455 list hence inherit the restriction used for neighbors. 

456 """ 

457 conn = get_neighbor_list(atoms) 

458 

459 if conn is None: 

460 conn = get_neighborlist(atoms) 

461 

462 no_of_loops = [0] * 8 

463 for s1 in range(len(atoms)): 

464 for s2 in conn[s1]: 

465 v12 = [s1] + [s2] 

466 for s3 in [s for s in conn[s2] if s not in v12]: 

467 v13 = v12 + [s3] 

468 if s1 in conn[s3]: 

469 no_of_loops[3] += 1 

470 for s4 in [s for s in conn[s3] if s not in v13]: 

471 v14 = v13 + [s4] 

472 if s1 in conn[s4]: 

473 no_of_loops[4] += 1 

474 for s5 in [s for s in conn[s4] if s not in v14]: 

475 v15 = v14 + [s5] 

476 if s1 in conn[s5]: 

477 no_of_loops[5] += 1 

478 for s6 in [s for s in conn[s5] if s not in v15]: 

479 v16 = v15 + [s6] 

480 if s1 in conn[s6]: 

481 no_of_loops[6] += 1 

482 for s7 in [s for s in conn[s6] if s not in v16]: 

483 # v17 = v16 + [s7] 

484 if s1 in conn[s7]: 

485 no_of_loops[7] += 1 

486 

487 to_return = [] 

488 for ring in rings: 

489 to_return.append(no_of_loops[ring]) 

490 

491 return to_return 

492 

493 

494def get_cell_angles_lengths(cell): 

495 """Returns cell vectors lengths (a,b,c) as well as different 

496 angles (alpha, beta, gamma, phi, chi, psi) (in radians). 

497 """ 

498 cellpar = cell_to_cellpar(cell) 

499 cellpar[3:] *= np.pi / 180 # convert angles to radians 

500 parnames = ['a', 'b', 'c', 'alpha', 'beta', 'gamma'] 

501 values = {n: p for n, p in zip(parnames, cellpar)} 

502 

503 volume = abs(np.linalg.det(cell)) 

504 for i, param in enumerate(['phi', 'chi', 'psi']): 

505 ab = np.linalg.norm( 

506 np.cross(cell[(i + 1) % 3, :], cell[(i + 2) % 3, :])) 

507 c = np.linalg.norm(cell[i, :]) 

508 s = np.abs(volume / (ab * c)) 

509 if 1 + 1e-6 > s > 1: 

510 s = 1. 

511 values[param] = np.arcsin(s) 

512 

513 return values 

514 

515 

516class CellBounds: 

517 """Class for defining as well as checking limits on 

518 cell vector lengths and angles. 

519 

520 Parameters: 

521 

522 bounds: dict 

523 Any of the following keywords can be used, in 

524 conjunction with a [low, high] list determining 

525 the lower and upper bounds: 

526 

527 a, b, c: 

528 Minimal and maximal lengths (in Angstrom) 

529 for the 1st, 2nd and 3rd lattice vectors. 

530 

531 alpha, beta, gamma: 

532 Minimal and maximal values (in degrees) 

533 for the angles between the lattice vectors. 

534 

535 phi, chi, psi: 

536 Minimal and maximal values (in degrees) 

537 for the angles between each lattice vector 

538 and the plane defined by the other two vectors. 

539 

540 Example: 

541 

542 >>> from ase.ga.utilities import CellBounds 

543 >>> CellBounds(bounds={'phi': [20, 160], 

544 ... 'chi': [60, 120], 

545 ... 'psi': [20, 160], 

546 ... 'a': [2, 20], 'b': [2, 20], 'c': [2, 20]}) 

547 """ 

548 

549 def __init__(self, bounds={}): 

550 self.bounds = {'alpha': [0, np.pi], 'beta': [0, np.pi], 

551 'gamma': [0, np.pi], 'phi': [0, np.pi], 

552 'chi': [0, np.pi], 'psi': [0, np.pi], 

553 'a': [0, 1e6], 'b': [0, 1e6], 'c': [0, 1e6]} 

554 

555 for param, bound in bounds.items(): 

556 if param not in ['a', 'b', 'c']: 

557 # Convert angle from degree to radians 

558 bound = [b * np.pi / 180. for b in bound] 

559 self.bounds[param] = bound 

560 

561 def is_within_bounds(self, cell): 

562 values = get_cell_angles_lengths(cell) 

563 verdict = True 

564 for param, bound in self.bounds.items(): 

565 if not (bound[0] <= values[param] <= bound[1]): 

566 verdict = False 

567 return verdict 

568 

569 

570def get_rotation_matrix(u, t): 

571 """Returns the transformation matrix for rotation over 

572 an angle t along an axis with direction u. 

573 """ 

574 ux, uy, uz = u 

575 cost, sint = np.cos(t), np.sin(t) 

576 rotmat = np.array([[(ux**2) * (1 - cost) + cost, 

577 ux * uy * (1 - cost) - uz * sint, 

578 ux * uz * (1 - cost) + uy * sint], 

579 [ux * uy * (1 - cost) + uz * sint, 

580 (uy**2) * (1 - cost) + cost, 

581 uy * uz * (1 - cost) - ux * sint], 

582 [ux * uz * (1 - cost) - uy * sint, 

583 uy * uz * (1 - cost) + ux * sint, 

584 (uz**2) * (1 - cost) + cost]]) 

585 return rotmat