Coverage for /builds/ase/ase/ase/ga/utilities.py: 55.24%
315 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""Various utility methods used troughout the GA."""
4import itertools
5import os
6import time
8import numpy as np
9from scipy.spatial.distance import cdist
11from ase.data import covalent_radii
12from ase.ga import get_neighbor_list
13from ase.geometry.cell import cell_to_cellpar
14from ase.geometry.rdf import get_rdf
15from ase.io import read, write
18def closest_distances_generator(atom_numbers, ratio_of_covalent_radii):
19 """Generates the blmin dict used across the GA.
20 The distances are based on the covalent radii of the atoms.
21 """
22 cr = covalent_radii
23 ratio = ratio_of_covalent_radii
25 blmin = {}
26 for i in atom_numbers:
27 blmin[(i, i)] = cr[i] * 2 * ratio
28 for j in atom_numbers:
29 if i == j:
30 continue
31 if (i, j) in blmin:
32 continue
33 blmin[(i, j)] = blmin[(j, i)] = ratio * (cr[i] + cr[j])
34 return blmin
37def get_mic_distance(p1, p2, cell, pbc):
38 """This method calculates the shortest distance between p1 and p2
39 through the cell boundaries defined by cell and pbc.
40 This method works for reasonable unit cells, but not for extremely
41 elongated ones.
42 """
43 ct = cell.T
44 pos = np.array((p1, p2))
45 scaled = np.linalg.solve(ct, pos.T).T
46 for i in range(3):
47 if pbc[i]:
48 scaled[:, i] %= 1.0
49 scaled[:, i] %= 1.0
50 P = np.dot(scaled, cell)
52 pbc_directions = [[-1, 1] * int(direction) + [0] for direction in pbc]
53 translations = np.array(list(itertools.product(*pbc_directions))).T
54 p0r = np.tile(np.reshape(P[0, :], (3, 1)), (1, translations.shape[1]))
55 p1r = np.tile(np.reshape(P[1, :], (3, 1)), (1, translations.shape[1]))
56 dp_vec = p0r + np.dot(ct, translations)
57 d = np.min(np.power(p1r - dp_vec, 2).sum(axis=0))**0.5
58 return d
61def db_call_with_error_tol(db_cursor, expression, args=[]):
62 """In case the GA is used on older versions of networking
63 filesystems there might be some delays. For this reason
64 some extra error tolerance when calling the SQLite db is
65 employed.
66 """
67 import sqlite3
68 for _ in range(10):
69 try:
70 db_cursor.execute(expression, args)
71 return
72 except sqlite3.OperationalError as e:
73 print(e)
74 time.sleep(2.)
75 raise sqlite3.OperationalError(
76 'Database still locked after 10 attempts (20 s)')
79def save_trajectory(confid, trajectory, folder):
80 """Saves traj files to the database folder.
81 This method should never be used directly,
82 but only through the DataConnection object.
83 """
84 fname = os.path.join(folder, 'traj%05d.traj' % confid)
85 write(fname, trajectory)
86 return fname
89def get_trajectory(fname):
90 """Extra error tolerance when loading traj files."""
91 fname = str(fname)
92 try:
93 t = read(fname)
94 except OSError as e:
95 print('get_trajectory error ' + e)
96 return t
99def gather_atoms_by_tag(atoms):
100 """Translates same-tag atoms so that they lie 'together',
101 with distance vectors as in the minimum image convention."""
102 tags = atoms.get_tags()
103 pos = atoms.get_positions()
104 for tag in list(set(tags)):
105 indices = np.where(tags == tag)[0]
106 if len(indices) == 1:
107 continue
108 vectors = atoms.get_distances(indices[0], indices[1:],
109 mic=True, vector=True)
110 pos[indices[1:]] = pos[indices[0]] + vectors
111 atoms.set_positions(pos)
114def atoms_too_close(atoms, bl, use_tags=False):
115 """Checks if any atoms in a are too close, as defined by
116 the distances in the bl dictionary.
118 use_tags: whether to use the Atoms tags to disable distance
119 checking within a set of atoms with the same tag.
121 Note: if certain atoms are constrained and use_tags is True,
122 this method may return unexpected results in case the
123 contraints prevent same-tag atoms to be gathered together in
124 the minimum-image-convention. In such cases, one should
125 (1) release the relevant constraints,
126 (2) apply the gather_atoms_by_tag function, and
127 (3) re-apply the constraints, before using the
128 atoms_too_close function.
129 """
130 a = atoms.copy()
131 if use_tags:
132 gather_atoms_by_tag(a)
134 pbc = a.get_pbc()
135 cell = a.get_cell()
136 num = a.get_atomic_numbers()
137 pos = a.get_positions()
138 tags = a.get_tags()
139 unique_types = sorted(list(set(num)))
141 neighbours = []
142 for i in range(3):
143 if pbc[i]:
144 neighbours.append([-1, 0, 1])
145 else:
146 neighbours.append([0])
148 for nx, ny, nz in itertools.product(*neighbours):
149 displacement = np.dot(cell.T, np.array([nx, ny, nz]).T)
150 pos_new = pos + displacement
151 distances = cdist(pos, pos_new)
153 if nx == 0 and ny == 0 and nz == 0:
154 if use_tags and len(a) > 1:
155 x = np.array([tags]).T
156 distances += 1e2 * (cdist(x, x) == 0)
157 else:
158 distances += 1e2 * np.identity(len(a))
160 iterator = itertools.combinations_with_replacement(unique_types, 2)
161 for type1, type2 in iterator:
162 x1 = np.where(num == type1)
163 x2 = np.where(num == type2)
164 if np.min(distances[x1].T[x2]) < bl[(type1, type2)]:
165 return True
167 return False
170def atoms_too_close_two_sets(a, b, bl):
171 """Checks if any atoms in a are too close to an atom in b,
172 as defined by the bl dictionary."""
173 pbc_a = a.get_pbc()
174 pbc_b = b.get_pbc()
175 cell_a = a.get_cell()
176 cell_b = a.get_cell()
177 assert np.allclose(pbc_a, pbc_b), (pbc_a, pbc_b)
178 assert np.allclose(cell_a, cell_b), (cell_a, cell_b)
180 pos_a = a.get_positions()
181 pos_b = b.get_positions()
183 num_a = a.get_atomic_numbers()
184 num_b = b.get_atomic_numbers()
185 unique_types = sorted(set(list(num_a) + list(num_b)))
187 neighbours = []
188 for i in range(3):
189 neighbours.append([-1, 0, 1] if pbc_a[i] else [0])
191 for nx, ny, nz in itertools.product(*neighbours):
192 displacement = np.dot(cell_a.T, np.array([nx, ny, nz]).T)
193 pos_b_disp = pos_b + displacement
194 distances = cdist(pos_a, pos_b_disp)
196 for type1 in unique_types:
197 if type1 not in num_a:
198 continue
199 x1 = np.where(num_a == type1)
201 for type2 in unique_types:
202 if type2 not in num_b:
203 continue
204 x2 = np.where(num_b == type2)
205 if np.min(distances[x1].T[x2]) < bl[(type1, type2)]:
206 return True
207 return False
210def get_all_atom_types(slab, atom_numbers_to_optimize):
211 """Utility method used to extract all unique atom types
212 from the atoms object slab and the list of atomic numbers
213 atom_numbers_to_optimize.
214 """
215 from_slab = list(set(slab.numbers))
216 from_top = list(set(atom_numbers_to_optimize))
217 from_slab.extend(from_top)
218 return list(set(from_slab))
221def get_distance_matrix(atoms, self_distance=1000):
222 """NB: This function is way slower than atoms.get_all_distances()
224 Returns a numpy matrix with the distances between the atoms
225 in the supplied atoms object, with the indices of the matrix
226 corresponding to the indices in the atoms object.
228 The parameter self_distance will be put in the diagonal
229 elements ([i][i])
230 """
231 dm = np.zeros([len(atoms), len(atoms)])
232 for i in range(len(atoms)):
233 dm[i][i] = self_distance
234 for j in range(i + 1, len(atoms)):
235 rij = atoms.get_distance(i, j)
236 dm[i][j] = rij
237 dm[j][i] = rij
238 return dm
241def get_nndist(atoms, distance_matrix):
242 """Returns an estimate of the nearest neighbor bond distance
243 in the supplied atoms object given the supplied distance_matrix.
245 The estimate comes from the first peak in the radial distribution
246 function.
247 """
248 rmax = 10. # No bonds longer than 10 angstrom expected
249 nbins = 200
250 rdf, dists = get_rdf(atoms, rmax, nbins, distance_matrix)
251 return dists[np.argmax(rdf)]
254def get_nnmat(atoms, mic=False):
255 """Calculate the nearest neighbor matrix as specified in
256 S. Lysgaard et al., Top. Catal., 2014, 57 (1-4), pp 33-39
258 Returns an array of average numbers of nearest neighbors
259 the order is determined by self.elements.
260 Example: self.elements = ["Cu", "Ni"]
261 get_nnmat returns a single list [Cu-Cu bonds/N(Cu),
262 Cu-Ni bonds/N(Cu), Ni-Cu bonds/N(Ni), Ni-Ni bonds/N(Ni)]
263 where N(element) is the number of atoms of the type element
264 in the atoms object.
266 The distance matrix can be quite costly to calculate every
267 time nnmat is required (and disk intensive if saved), thus
268 it makes sense to calculate nnmat along with e.g. the
269 potential energy and save it in atoms.info['data']['nnmat'].
270 """
271 if 'data' in atoms.info and 'nnmat' in atoms.info['data']:
272 return atoms.info['data']['nnmat']
273 elements = sorted(set(atoms.get_chemical_symbols()))
274 nnmat = np.zeros((len(elements), len(elements)))
275 # dm = get_distance_matrix(atoms)
276 dm = atoms.get_all_distances(mic=mic)
277 nndist = get_nndist(atoms, dm) + 0.2
278 for i in range(len(atoms)):
279 row = [j for j in range(len(elements))
280 if atoms[i].symbol == elements[j]][0]
281 neighbors = [j for j in range(len(dm[i])) if dm[i][j] < nndist]
282 for n in neighbors:
283 column = [j for j in range(len(elements))
284 if atoms[n].symbol == elements[j]][0]
285 nnmat[row][column] += 1
286 # divide by the number of that type of atoms in the structure
287 for i, el in enumerate(elements):
288 nnmat[i] /= len([j for j in range(len(atoms))
289 if atoms[int(j)].symbol == el])
290 # makes a single list out of a list of lists
291 nnlist = np.reshape(nnmat, (len(nnmat)**2))
292 return nnlist
295def get_nnmat_string(atoms, decimals=2, mic=False):
296 nnmat = get_nnmat(atoms, mic=mic)
297 s = '-'.join(['{1:2.{0}f}'.format(decimals, i)
298 for i in nnmat])
299 if len(nnmat) == 1:
300 return s + '-'
301 return s
304def get_connections_index(atoms, max_conn=5, no_count_types=None):
305 """This method returns a dictionary where each key value are a
306 specific number of neighbors and list of atoms indices with
307 that amount of neighbors respectively. The method utilizes the
308 neighbor list and hence inherit the restrictions for
309 neighbors. Option added to remove connections between
310 defined atom types.
312 Parameters
313 ----------
315 atoms : Atoms object
316 The connections will be counted using this supplied Atoms object
318 max_conn : int
319 Any atom with more connections than this will be counted as
320 having max_conn connections.
321 Default 5
323 no_count_types : list or None
324 List of atomic numbers that should be excluded in the count.
325 Default None (meaning all atoms count).
326 """
327 conn = get_neighbor_list(atoms)
329 if conn is None:
330 conn = get_neighborlist(atoms)
332 if no_count_types is None:
333 no_count_types = []
335 conn_index = {}
336 for i in range(len(atoms)):
337 if atoms[i].number not in no_count_types:
338 cconn = min(len(conn[i]), max_conn - 1)
339 if cconn not in conn_index:
340 conn_index[cconn] = []
341 conn_index[cconn].append(i)
343 return conn_index
346def get_atoms_connections(atoms, max_conn=5, no_count_types=None):
347 """This method returns a list of the numbers of atoms
348 with X number of neighbors. The method utilizes the
349 neighbor list and hence inherit the restrictions for
350 neighbors. Option added to remove connections between
351 defined atom types.
352 """
353 conn_index = get_connections_index(atoms, max_conn=max_conn,
354 no_count_types=no_count_types)
356 no_of_conn = [0] * max_conn
357 for i in conn_index:
358 no_of_conn[i] += len(conn_index[i])
360 return no_of_conn
363def get_angles_distribution(atoms, ang_grid=9):
364 """Method to get the distribution of bond angles
365 in bins (default 9) with bonds defined from
366 the get_neighbor_list().
367 """
368 conn = get_neighbor_list(atoms)
370 if conn is None:
371 conn = get_neighborlist(atoms)
373 bins = [0] * ang_grid
375 for atom in atoms:
376 for i in conn[atom.index]:
377 for j in conn[atom.index]:
378 if j != i:
379 a = atoms.get_angle(i, atom.index, j)
380 for k in range(ang_grid):
381 if (k + 1) * 180. / ang_grid > a > k * 180. / ang_grid:
382 bins[k] += 1
383 # Removing dobbelt counting
384 for i in range(ang_grid):
385 bins[i] /= 2.
386 return bins
389def get_neighborlist(atoms, dx=0.2, no_count_types=None):
390 """Method to get the a dict with list of neighboring
391 atoms defined as the two covalent radii + fixed distance.
392 Option added to remove neighbors between defined atom types.
393 """
394 cell = atoms.get_cell()
395 pbc = atoms.get_pbc()
397 if no_count_types is None:
398 no_count_types = []
400 conn = {}
401 for atomi in atoms:
402 conn_this_atom = []
403 for atomj in atoms:
404 if atomi.index != atomj.index:
405 if atomi.number not in no_count_types:
406 if atomj.number not in no_count_types:
407 d = get_mic_distance(atomi.position,
408 atomj.position,
409 cell,
410 pbc)
411 cri = covalent_radii[atomi.number]
412 crj = covalent_radii[atomj.number]
413 d_max = crj + cri + dx
414 if d < d_max:
415 conn_this_atom.append(atomj.index)
416 conn[atomi.index] = conn_this_atom
417 return conn
420def get_atoms_distribution(atoms, number_of_bins=5, max_distance=8,
421 center=None, no_count_types=None):
422 """Method to get the distribution of atoms in the
423 structure in bins of distances from a defined
424 center. Option added to remove counting of
425 certain atom types.
426 """
427 pbc = atoms.get_pbc()
428 cell = atoms.get_cell()
429 if center is None:
430 # Center used for the atom distribution if None is supplied!
431 cx = sum(cell[:, 0]) / 2.
432 cy = sum(cell[:, 1]) / 2.
433 cz = sum(cell[:, 2]) / 2.
434 center = (cx, cy, cz)
435 bins = [0] * number_of_bins
437 if no_count_types is None:
438 no_count_types = []
440 for atom in atoms:
441 if atom.number not in no_count_types:
442 d = get_mic_distance(atom.position, center, cell, pbc)
443 for k in range(number_of_bins - 1):
444 min_dis_cur_bin = k * max_distance / (number_of_bins - 1.)
445 max_dis_cur_bin = ((k + 1) * max_distance /
446 (number_of_bins - 1.))
447 if min_dis_cur_bin < d < max_dis_cur_bin:
448 bins[k] += 1
449 if d > max_distance:
450 bins[number_of_bins - 1] += 1
451 return bins
454def get_rings(atoms, rings=[5, 6, 7]):
455 """This method return a list of the number of atoms involved
456 in rings in the structures. It uses the neighbor
457 list hence inherit the restriction used for neighbors.
458 """
459 conn = get_neighbor_list(atoms)
461 if conn is None:
462 conn = get_neighborlist(atoms)
464 no_of_loops = [0] * 8
465 for s1 in range(len(atoms)):
466 for s2 in conn[s1]:
467 v12 = [s1] + [s2]
468 for s3 in [s for s in conn[s2] if s not in v12]:
469 v13 = v12 + [s3]
470 if s1 in conn[s3]:
471 no_of_loops[3] += 1
472 for s4 in [s for s in conn[s3] if s not in v13]:
473 v14 = v13 + [s4]
474 if s1 in conn[s4]:
475 no_of_loops[4] += 1
476 for s5 in [s for s in conn[s4] if s not in v14]:
477 v15 = v14 + [s5]
478 if s1 in conn[s5]:
479 no_of_loops[5] += 1
480 for s6 in [s for s in conn[s5] if s not in v15]:
481 v16 = v15 + [s6]
482 if s1 in conn[s6]:
483 no_of_loops[6] += 1
484 for s7 in [s for s in conn[s6] if s not in v16]:
485 # v17 = v16 + [s7]
486 if s1 in conn[s7]:
487 no_of_loops[7] += 1
489 to_return = []
490 for ring in rings:
491 to_return.append(no_of_loops[ring])
493 return to_return
496def get_cell_angles_lengths(cell):
497 """Returns cell vectors lengths (a,b,c) as well as different
498 angles (alpha, beta, gamma, phi, chi, psi) (in radians).
499 """
500 cellpar = cell_to_cellpar(cell)
501 cellpar[3:] *= np.pi / 180 # convert angles to radians
502 parnames = ['a', 'b', 'c', 'alpha', 'beta', 'gamma']
503 values = {n: p for n, p in zip(parnames, cellpar)}
505 volume = abs(np.linalg.det(cell))
506 for i, param in enumerate(['phi', 'chi', 'psi']):
507 ab = np.linalg.norm(
508 np.cross(cell[(i + 1) % 3, :], cell[(i + 2) % 3, :]))
509 c = np.linalg.norm(cell[i, :])
510 s = np.abs(volume / (ab * c))
511 if 1 + 1e-6 > s > 1:
512 s = 1.
513 values[param] = np.arcsin(s)
515 return values
518class CellBounds:
519 """Class for defining as well as checking limits on
520 cell vector lengths and angles.
522 Parameters:
524 bounds: dict
525 Any of the following keywords can be used, in
526 conjunction with a [low, high] list determining
527 the lower and upper bounds:
529 a, b, c:
530 Minimal and maximal lengths (in Angstrom)
531 for the 1st, 2nd and 3rd lattice vectors.
533 alpha, beta, gamma:
534 Minimal and maximal values (in degrees)
535 for the angles between the lattice vectors.
537 phi, chi, psi:
538 Minimal and maximal values (in degrees)
539 for the angles between each lattice vector
540 and the plane defined by the other two vectors.
542 Example:
544 >>> from ase.ga.utilities import CellBounds
545 >>> CellBounds(bounds={'phi': [20, 160],
546 ... 'chi': [60, 120],
547 ... 'psi': [20, 160],
548 ... 'a': [2, 20], 'b': [2, 20], 'c': [2, 20]})
549 """
551 def __init__(self, bounds={}):
552 self.bounds = {'alpha': [0, np.pi], 'beta': [0, np.pi],
553 'gamma': [0, np.pi], 'phi': [0, np.pi],
554 'chi': [0, np.pi], 'psi': [0, np.pi],
555 'a': [0, 1e6], 'b': [0, 1e6], 'c': [0, 1e6]}
557 for param, bound in bounds.items():
558 if param not in ['a', 'b', 'c']:
559 # Convert angle from degree to radians
560 bound = [b * np.pi / 180. for b in bound]
561 self.bounds[param] = bound
563 def is_within_bounds(self, cell):
564 values = get_cell_angles_lengths(cell)
565 verdict = True
566 for param, bound in self.bounds.items():
567 if not (bound[0] <= values[param] <= bound[1]):
568 verdict = False
569 return verdict
572def get_rotation_matrix(u, t):
573 """Returns the transformation matrix for rotation over
574 an angle t along an axis with direction u.
575 """
576 ux, uy, uz = u
577 cost, sint = np.cos(t), np.sin(t)
578 rotmat = np.array([[(ux**2) * (1 - cost) + cost,
579 ux * uy * (1 - cost) - uz * sint,
580 ux * uz * (1 - cost) + uy * sint],
581 [ux * uy * (1 - cost) + uz * sint,
582 (uy**2) * (1 - cost) + cost,
583 uy * uz * (1 - cost) - ux * sint],
584 [ux * uz * (1 - cost) - uy * sint,
585 uy * uz * (1 - cost) + ux * sint,
586 (uz**2) * (1 - cost) + cost]])
587 return rotmat