Coverage for /builds/debichem-team/python-ase/ase/ga/utilities.py : 57.70%
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 math
5import itertools
6import numpy as np
7from scipy.spatial.distance import cdist
8from ase.io import write, read
9from ase.geometry.cell import cell_to_cellpar
10from ase.data import covalent_radii
11from ase.ga import get_neighbor_list
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
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
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)
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
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)')
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 IOError 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)
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.
116 use_tags: whether to use the Atoms tags to disable distance
117 checking within a set of atoms with the same tag.
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)
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)))
139 neighbours = []
140 for i in range(3):
141 if pbc[i]:
142 neighbours.append([-1, 0, 1])
143 else:
144 neighbours.append([0])
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)
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))
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
165 return False
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)
178 pos_a = a.get_positions()
179 pos_b = b.get_positions()
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)))
185 neighbours = []
186 for i in range(3):
187 neighbours.append([-1, 0, 1] if pbc_a[i] else [0])
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)
194 for type1 in unique_types:
195 if type1 not in num_a:
196 continue
197 x1 = np.where(num_a == type1)
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
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))
219def get_distance_matrix(atoms, self_distance=1000):
220 """NB: This function is way slower than atoms.get_all_distances()
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.
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
239def get_rdf(atoms, rmax, nbins, distance_matrix=None,
240 elements=None, no_dists=False):
241 """Returns two numpy arrays; the radial distribution function
242 and the corresponding distances of the supplied atoms object.
243 If no_dists = True then only the first array is returned.
245 Note that the rdf is computed following the standard solid state
246 definition which uses the cell volume in the normalization.
247 This may or may not be appropriate in cases where one or more
248 directions is non-periodic.
250 Parameters:
252 rmax : float
253 The maximum distance that will contribute to the rdf.
254 The unit cell should be large enough so that it encloses a
255 sphere with radius rmax in the periodic directions.
257 nbins : int
258 Number of bins to divide the rdf into.
260 distance_matrix : numpy.array
261 An array of distances between atoms, typically
262 obtained by atoms.get_all_distances().
263 Default None meaning that it will be calculated.
265 elements : list or tuple
266 List of two atomic numbers. If elements is not None the partial
267 rdf for the supplied elements will be returned.
269 no_dists : bool
270 If True then the second array with rdf distances will not be returned
271 """
272 # First check whether the cell is sufficiently large
273 cell = atoms.get_cell()
274 vol = atoms.get_volume()
275 pbc = atoms.get_pbc()
276 for i in range(3):
277 if pbc[i]:
278 axb = np.cross(cell[(i + 1) % 3, :], cell[(i + 2) % 3, :])
279 h = vol / np.linalg.norm(axb)
280 assert h > 2 * rmax, 'The cell is not large enough in ' \
281 'direction %d: %.3f < 2*rmax=%.3f' % (i, h, 2 * rmax)
283 dm = distance_matrix
284 if dm is None:
285 dm = atoms.get_all_distances(mic=True)
286 rdf = np.zeros(nbins + 1)
287 dr = float(rmax / nbins)
289 if elements is None:
290 # Coefficients to use for normalization
291 phi = len(atoms) / vol
292 norm = 2.0 * math.pi * dr * phi * len(atoms)
294 for i in range(len(atoms)):
295 for j in range(i + 1, len(atoms)):
296 rij = dm[i][j]
297 index = int(math.ceil(rij / dr))
298 if index <= nbins:
299 rdf[index] += 1
300 else:
301 i_indices = np.where(atoms.numbers == elements[0])[0]
302 phi = len(i_indices) / vol
303 norm = 4.0 * math.pi * dr * phi * len(atoms)
305 for i in i_indices:
306 for j in np.where(atoms.numbers == elements[1])[0]:
307 rij = dm[i][j]
308 index = int(math.ceil(rij / dr))
309 if index <= nbins:
310 rdf[index] += 1
312 dists = []
313 for i in range(1, nbins + 1):
314 rrr = (i - 0.5) * dr
315 dists.append(rrr)
316 # Normalize
317 rdf[i] /= (norm * ((rrr**2) + (dr**2) / 12.))
319 if no_dists:
320 return rdf[1:]
321 return rdf[1:], np.array(dists)
324def get_nndist(atoms, distance_matrix):
325 """Returns an estimate of the nearest neighbor bond distance
326 in the supplied atoms object given the supplied distance_matrix.
328 The estimate comes from the first peak in the radial distribution
329 function.
330 """
331 rmax = 10. # No bonds longer than 10 angstrom expected
332 nbins = 200
333 rdf, dists = get_rdf(atoms, rmax, nbins, distance_matrix)
334 return dists[np.argmax(rdf)]
337def get_nnmat(atoms, mic=False):
338 """Calculate the nearest neighbor matrix as specified in
339 S. Lysgaard et al., Top. Catal., 2014, 57 (1-4), pp 33-39
341 Returns an array of average numbers of nearest neighbors
342 the order is determined by self.elements.
343 Example: self.elements = ["Cu", "Ni"]
344 get_nnmat returns a single list [Cu-Cu bonds/N(Cu),
345 Cu-Ni bonds/N(Cu), Ni-Cu bonds/N(Ni), Ni-Ni bonds/N(Ni)]
346 where N(element) is the number of atoms of the type element
347 in the atoms object.
349 The distance matrix can be quite costly to calculate every
350 time nnmat is required (and disk intensive if saved), thus
351 it makes sense to calculate nnmat along with e.g. the
352 potential energy and save it in atoms.info['data']['nnmat'].
353 """
354 if 'data' in atoms.info and 'nnmat' in atoms.info['data']:
355 return atoms.info['data']['nnmat']
356 elements = sorted(set(atoms.get_chemical_symbols()))
357 nnmat = np.zeros((len(elements), len(elements)))
358 # dm = get_distance_matrix(atoms)
359 dm = atoms.get_all_distances(mic=mic)
360 nndist = get_nndist(atoms, dm) + 0.2
361 for i in range(len(atoms)):
362 row = [j for j in range(len(elements))
363 if atoms[i].symbol == elements[j]][0]
364 neighbors = [j for j in range(len(dm[i])) if dm[i][j] < nndist]
365 for n in neighbors:
366 column = [j for j in range(len(elements))
367 if atoms[n].symbol == elements[j]][0]
368 nnmat[row][column] += 1
369 # divide by the number of that type of atoms in the structure
370 for i, el in enumerate(elements):
371 nnmat[i] /= len([j for j in range(len(atoms))
372 if atoms[int(j)].symbol == el])
373 # makes a single list out of a list of lists
374 nnlist = np.reshape(nnmat, (len(nnmat)**2))
375 return nnlist
378def get_nnmat_string(atoms, decimals=2, mic=False):
379 nnmat = get_nnmat(atoms, mic=mic)
380 s = '-'.join(['{1:2.{0}f}'.format(decimals, i)
381 for i in nnmat])
382 if len(nnmat) == 1:
383 return s + '-'
384 return s
387def get_connections_index(atoms, max_conn=5, no_count_types=None):
388 """This method returns a dictionary where each key value are a
389 specific number of neighbors and list of atoms indices with
390 that amount of neighbors respectively. The method utilizes the
391 neighbor list and hence inherit the restrictions for
392 neighbors. Option added to remove connections between
393 defined atom types.
395 Parameters
396 ----------
398 atoms : Atoms object
399 The connections will be counted using this supplied Atoms object
401 max_conn : int
402 Any atom with more connections than this will be counted as
403 having max_conn connections.
404 Default 5
406 no_count_types : list or None
407 List of atomic numbers that should be excluded in the count.
408 Default None (meaning all atoms count).
409 """
410 conn = get_neighbor_list(atoms)
412 if conn is None:
413 conn = get_neighborlist(atoms)
415 if no_count_types is None:
416 no_count_types = []
418 conn_index = {}
419 for i in range(len(atoms)):
420 if atoms[i].number not in no_count_types:
421 cconn = min(len(conn[i]), max_conn - 1)
422 if cconn not in conn_index:
423 conn_index[cconn] = []
424 conn_index[cconn].append(i)
426 return conn_index
429def get_atoms_connections(atoms, max_conn=5, no_count_types=None):
430 """This method returns a list of the numbers of atoms
431 with X number of neighbors. The method utilizes the
432 neighbor list and hence inherit the restrictions for
433 neighbors. Option added to remove connections between
434 defined atom types.
435 """
436 conn_index = get_connections_index(atoms, max_conn=max_conn,
437 no_count_types=no_count_types)
439 no_of_conn = [0] * max_conn
440 for i in conn_index:
441 no_of_conn[i] += len(conn_index[i])
443 return no_of_conn
446def get_angles_distribution(atoms, ang_grid=9):
447 """Method to get the distribution of bond angles
448 in bins (default 9) with bonds defined from
449 the get_neighbor_list().
450 """
451 conn = get_neighbor_list(atoms)
453 if conn is None:
454 conn = get_neighborlist(atoms)
456 bins = [0] * ang_grid
458 for atom in atoms:
459 for i in conn[atom.index]:
460 for j in conn[atom.index]:
461 if j != i:
462 a = atoms.get_angle(i, atom.index, j)
463 for k in range(ang_grid):
464 if (k + 1) * 180. / ang_grid > a > k * 180. / ang_grid:
465 bins[k] += 1
466 # Removing dobbelt counting
467 for i in range(ang_grid):
468 bins[i] /= 2.
469 return bins
472def get_neighborlist(atoms, dx=0.2, no_count_types=None):
473 """Method to get the a dict with list of neighboring
474 atoms defined as the two covalent radii + fixed distance.
475 Option added to remove neighbors between defined atom types.
476 """
477 cell = atoms.get_cell()
478 pbc = atoms.get_pbc()
480 if no_count_types is None:
481 no_count_types = []
483 conn = {}
484 for atomi in atoms:
485 conn_this_atom = []
486 for atomj in atoms:
487 if atomi.index != atomj.index:
488 if atomi.number not in no_count_types:
489 if atomj.number not in no_count_types:
490 d = get_mic_distance(atomi.position,
491 atomj.position,
492 cell,
493 pbc)
494 cri = covalent_radii[atomi.number]
495 crj = covalent_radii[atomj.number]
496 d_max = crj + cri + dx
497 if d < d_max:
498 conn_this_atom.append(atomj.index)
499 conn[atomi.index] = conn_this_atom
500 return conn
503def get_atoms_distribution(atoms, number_of_bins=5, max_distance=8,
504 center=None, no_count_types=None):
505 """Method to get the distribution of atoms in the
506 structure in bins of distances from a defined
507 center. Option added to remove counting of
508 certain atom types.
509 """
510 pbc = atoms.get_pbc()
511 cell = atoms.get_cell()
512 if center is None:
513 # Center used for the atom distribution if None is supplied!
514 cx = sum(cell[:, 0]) / 2.
515 cy = sum(cell[:, 1]) / 2.
516 cz = sum(cell[:, 2]) / 2.
517 center = (cx, cy, cz)
518 bins = [0] * number_of_bins
520 if no_count_types is None:
521 no_count_types = []
523 for atom in atoms:
524 if atom.number not in no_count_types:
525 d = get_mic_distance(atom.position, center, cell, pbc)
526 for k in range(number_of_bins - 1):
527 min_dis_cur_bin = k * max_distance / (number_of_bins - 1.)
528 max_dis_cur_bin = ((k + 1) * max_distance /
529 (number_of_bins - 1.))
530 if min_dis_cur_bin < d < max_dis_cur_bin:
531 bins[k] += 1
532 if d > max_distance:
533 bins[number_of_bins - 1] += 1
534 return bins
537def get_rings(atoms, rings=[5, 6, 7]):
538 """This method return a list of the number of atoms involved
539 in rings in the structures. It uses the neighbor
540 list hence inherit the restriction used for neighbors.
541 """
542 conn = get_neighbor_list(atoms)
544 if conn is None:
545 conn = get_neighborlist(atoms)
547 no_of_loops = [0] * 8
548 for s1 in range(len(atoms)):
549 for s2 in conn[s1]:
550 v12 = [s1] + [s2]
551 for s3 in [s for s in conn[s2] if s not in v12]:
552 v13 = v12 + [s3]
553 if s1 in conn[s3]:
554 no_of_loops[3] += 1
555 for s4 in [s for s in conn[s3] if s not in v13]:
556 v14 = v13 + [s4]
557 if s1 in conn[s4]:
558 no_of_loops[4] += 1
559 for s5 in [s for s in conn[s4] if s not in v14]:
560 v15 = v14 + [s5]
561 if s1 in conn[s5]:
562 no_of_loops[5] += 1
563 for s6 in [s for s in conn[s5] if s not in v15]:
564 v16 = v15 + [s6]
565 if s1 in conn[s6]:
566 no_of_loops[6] += 1
567 for s7 in [s for s in conn[s6] if s not in v16]:
568 # v17 = v16 + [s7]
569 if s1 in conn[s7]:
570 no_of_loops[7] += 1
572 to_return = []
573 for ring in rings:
574 to_return.append(no_of_loops[ring])
576 return to_return
579def get_cell_angles_lengths(cell):
580 """Returns cell vectors lengths (a,b,c) as well as different
581 angles (alpha, beta, gamma, phi, chi, psi) (in radians).
582 """
583 cellpar = cell_to_cellpar(cell)
584 cellpar[3:] *= np.pi / 180 # convert angles to radians
585 parnames = ['a', 'b', 'c', 'alpha', 'beta', 'gamma']
586 values = {n: p for n, p in zip(parnames, cellpar)}
588 volume = abs(np.linalg.det(cell))
589 for i, param in enumerate(['phi', 'chi', 'psi']):
590 ab = np.linalg.norm(
591 np.cross(cell[(i + 1) % 3, :], cell[(i + 2) % 3, :]))
592 c = np.linalg.norm(cell[i, :])
593 s = np.abs(volume / (ab * c))
594 if 1 + 1e-6 > s > 1:
595 s = 1.
596 values[param] = np.arcsin(s)
598 return values
601class CellBounds:
602 """Class for defining as well as checking limits on
603 cell vector lengths and angles.
605 Parameters:
607 bounds: dict
608 Any of the following keywords can be used, in
609 conjunction with a [low, high] list determining
610 the lower and upper bounds:
612 a, b, c:
613 Minimal and maximal lengths (in Angstrom)
614 for the 1st, 2nd and 3rd lattice vectors.
616 alpha, beta, gamma:
617 Minimal and maximal values (in degrees)
618 for the angles between the lattice vectors.
620 phi, chi, psi:
621 Minimal and maximal values (in degrees)
622 for the angles between each lattice vector
623 and the plane defined by the other two vectors.
625 Example:
627 >>> from ase.ga.utilities import CellBounds
628 >>> CellBounds(bounds={'phi': [20, 160],
629 ... 'chi': [60, 120],
630 ... 'psi': [20, 160],
631 ... 'a': [2, 20], 'b': [2, 20], 'c': [2, 20]})
632 """
633 def __init__(self, bounds={}):
634 self.bounds = {'alpha': [0, np.pi], 'beta': [0, np.pi],
635 'gamma': [0, np.pi], 'phi': [0, np.pi],
636 'chi': [0, np.pi], 'psi': [0, np.pi],
637 'a': [0, 1e6], 'b': [0, 1e6], 'c': [0, 1e6]}
639 for param, bound in bounds.items():
640 if param not in ['a', 'b', 'c']:
641 # Convert angle from degree to radians
642 bound = [b * np.pi / 180. for b in bound]
643 self.bounds[param] = bound
645 def is_within_bounds(self, cell):
646 values = get_cell_angles_lengths(cell)
647 verdict = True
648 for param, bound in self.bounds.items():
649 if not (bound[0] <= values[param] <= bound[1]):
650 verdict = False
651 return verdict
654def get_rotation_matrix(u, t):
655 """Returns the transformation matrix for rotation over
656 an angle t along an axis with direction u.
657 """
658 ux, uy, uz = u
659 cost, sint = np.cos(t), np.sin(t)
660 rotmat = np.array([[(ux**2) * (1 - cost) + cost,
661 ux * uy * (1 - cost) - uz * sint,
662 ux * uz * (1 - cost) + uy * sint],
663 [ux * uy * (1 - cost) + uz * sint,
664 (uy**2) * (1 - cost) + cost,
665 uy * uz * (1 - cost) - ux * sint],
666 [ux * uz * (1 - cost) - uy * sint,
667 uy * uz * (1 - cost) + ux * sint,
668 (uz**2) * (1 - cost) + cost]])
669 return rotmat