Coverage for /builds/kinetik161/ase/ase/ga/utilities.py: 52.70%

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

2import itertools

3import os

4import time

6import numpy as np

7from scipy.spatial.distance import cdist

9from ase.data import covalent_radii

10from ase.ga import get_neighbor_list

11from ase.geometry.cell import cell_to_cellpar

12from ase.geometry.rdf import get_rdf

13from ase.io import read, write

16def closest_distances_generator(atom_numbers, ratio_of_covalent_radii):

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

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

19 """

20 cr = covalent_radii

21 ratio = ratio_of_covalent_radii

23 blmin = {}

24 for i in atom_numbers:

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

26 for j in atom_numbers:

27 if i == j:

28 continue

29 if (i, j) in blmin:

30 continue

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

32 return blmin

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

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

37 through the cell boundaries defined by cell and pbc.

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

39 elongated ones.

40 """

41 ct = cell.T

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

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

44 for i in range(3):

45 if pbc[i]:

46 scaled[:, i] %= 1.0

47 scaled[:, i] %= 1.0

48 P = np.dot(scaled, cell)

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

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

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

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

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

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

56 return d

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

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

61 filesystems there might be some delays. For this reason

62 some extra error tolerance when calling the SQLite db is

63 employed.

64 """

65 import sqlite3

66 for i in range(10):

67 try:

68 db_cursor.execute(expression, args)

69 return

70 except sqlite3.OperationalError as e:

71 print(e)

72 time.sleep(2.)

73 raise sqlite3.OperationalError(

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

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

87def get_trajectory(fname):

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

89 fname = str(fname)

90 try:

91 t = read(fname)

92 except OSError as e:

93 print('get_trajectory error ' + e)

94 return t

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