Coverage for /builds/debichem-team/python-ase/ase/geometry/analysis.py : 63.75%
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# flake8: noqa
2"""Tools for analyzing instances of :class:`~ase.Atoms`
3"""
5from ase.neighborlist import build_neighbor_list, get_distance_matrix, get_distance_indices
6from ase.ga.utilities import get_rdf
7from ase import Atoms
10__all__ = ['Analysis']
13class Analysis:
14 """Analysis class
16 Parameters for initialization:
18 images: :class:`~ase.Atoms` object or list of such
19 Images to analyze.
20 nl: None, :class:`~ase.neighborlist.NeighborList` object or list of such
21 Neighborlist(s) for the given images. One or nImages, depending if bonding
22 pattern changes or is constant. Using one Neigborlist greatly improves speed.
23 kwargs: options, dict
24 Arguments for constructing :class:`~ase.neighborlist.NeighborList` object if :data:`nl` is None.
26 The choice of ``bothways=True`` for the :class:`~ase.neighborlist.NeighborList` object
27 will not influence the amount of bonds/angles/dihedrals you get, all are reported
28 in both directions. Use the *unique*-labeled properties to get lists without
29 duplicates.
30 """
32 def __init__(self, images, nl=None, **kwargs):
33 self.images = images
35 if isinstance(nl, list):
36 assert len(nl) == self.nImages
37 self._nl = nl
38 elif nl is not None:
39 self._nl = [ nl ]
40 else:
41 self._nl = [ build_neighbor_list(self.images[0], **kwargs) ]
43 self._cache = {}
45 def _get_slice(self, imageIdx):
46 """Return a slice from user input.
47 Using *imageIdx* (can be integer or slice) the analyzed frames can be specified.
48 If *imageIdx* is None, all frames will be analyzed.
49 """
50 #get slice from imageIdx
51 if isinstance(imageIdx, int):
52 sl = slice(imageIdx, imageIdx+1)
53 elif isinstance(imageIdx, slice):
54 sl = imageIdx
55 elif imageIdx is None:
56 sl = slice(0, None)
57 else:
58 raise ValueError("Unsupported type for imageIdx in ase.geometry.analysis.Analysis._get_slice")
59 return sl
61 @property
62 def images(self):
63 """Images.
65 Set during initialization but can also be set later.
66 """
67 return self._images
69 @images.setter
70 def images(self, images):
71 """Set images"""
72 if isinstance(images, list):
73 self._images = images
74 else:
75 self._images = [ images ]
78 @images.deleter
79 def images(self):
80 """Delete images"""
81 self._images = None
83 @property
84 def nImages(self):
85 """Number of Images in this instance.
87 Cannot be set, is determined automatically.
88 """
89 return len(self.images)
91 @property
92 def nl(self):
93 """Neighbor Lists in this instance.
95 Set during initialization.
97 **No setter or deleter, only getter**
98 """
99 return self._nl
101 def _get_all_x(self, distance):
102 """Helper function to get bonds, angles, dihedrals"""
103 maxIter = self.nImages
104 if len(self.nl) == 1:
105 maxIter = 1
107 xList = []
108 for i in range(maxIter):
109 xList.append(get_distance_indices(self.distance_matrix[i], distance))
111 return xList
113 @property
114 def all_bonds(self):
115 """All Bonds.
117 A list with indices of bonded atoms for each neighborlist in *self*.
118 Atom i is connected to all atoms inside result[i]. Duplicates from PBCs are
119 removed. See also :data:`unique_bonds`.
121 **No setter or deleter, only getter**
122 """
123 if not 'allBonds' in self._cache:
124 self._cache['allBonds'] = self._get_all_x(1)
126 return self._cache['allBonds']
128 @property
129 def all_angles(self):
130 """All angles
132 A list with indices of atoms in angles for each neighborlist in *self*.
133 Atom i forms an angle to the atoms inside the tuples in result[i]:
134 i -- result[i][x][0] -- result[i][x][1]
135 where x is in range(number of angles from i). See also :data:`unique_angles`.
137 **No setter or deleter, only getter**
138 """
139 if not 'allAngles' in self._cache:
140 self._cache['allAngles'] = []
141 distList = self._get_all_x(2)
143 for imI in range(len(distList)):
144 self._cache['allAngles'].append([])
145 #iterate over second neighbors of all atoms
146 for iAtom, secNeighs in enumerate(distList[imI]):
147 self._cache['allAngles'][-1].append([])
148 if len(secNeighs) == 0:
149 continue
150 firstNeighs = self.all_bonds[imI][iAtom]
151 #iterate over second neighbors of iAtom
152 for kAtom in secNeighs:
153 relevantFirstNeighs = [ idx for idx in firstNeighs if kAtom in self.all_bonds[imI][idx] ]
154 #iterate over all atoms that are connected to iAtom and kAtom
155 for jAtom in relevantFirstNeighs:
156 self._cache['allAngles'][-1][-1].append((jAtom, kAtom))
158 return self._cache['allAngles']
160 @property
161 def all_dihedrals(self):
162 """All dihedrals
164 Returns a list with indices of atoms in dihedrals for each neighborlist in this instance.
165 Atom i forms a dihedral to the atoms inside the tuples in result[i]:
166 i -- result[i][x][0] -- result[i][x][1] -- result[i][x][2]
167 where x is in range(number of dihedrals from i). See also :data:`unique_dihedrals`.
169 **No setter or deleter, only getter**
170 """
171 if not 'allDihedrals' in self._cache:
172 self._cache['allDihedrals'] = []
173 distList = self._get_all_x(3)
175 for imI in range(len(distList)):
176 self._cache['allDihedrals'].append([])
177 for iAtom, thirdNeighs in enumerate(distList[imI]):
178 self._cache['allDihedrals'][-1].append([])
179 if len(thirdNeighs) == 0:
180 continue
181 anglesI = self.all_angles[imI][iAtom]
182 #iterate over third neighbors of iAtom
183 for lAtom in thirdNeighs:
184 secondNeighs = [ angle[-1] for angle in anglesI ]
185 firstNeighs = [ angle[0] for angle in anglesI ]
186 relevantSecondNeighs = [ idx for idx in secondNeighs if lAtom in self.all_bonds[imI][idx] ]
187 relevantFirstNeighs = [ firstNeighs[secondNeighs.index(idx)] for idx in relevantSecondNeighs ]
188 #iterate over all atoms that are connected to iAtom and lAtom
189 for jAtom, kAtom in zip(relevantFirstNeighs, relevantSecondNeighs):
190 #remove dihedrals in circles
191 tupl = (jAtom, kAtom, lAtom)
192 if len(set((iAtom, ) + tupl)) != 4:
193 continue
194 #avoid duplicates
195 elif tupl in self._cache['allDihedrals'][-1][-1]:
196 continue
197 elif iAtom in tupl:
198 raise RuntimeError("Something is wrong in analysis.all_dihedrals!")
199 self._cache['allDihedrals'][-1][-1].append((jAtom, kAtom, lAtom))
201 return self._cache['allDihedrals']
203 @property
204 def adjacency_matrix(self):
205 """The adjacency/connectivity matrix.
207 If not already done, build a list of adjacency matrices for all :data:`nl`.
209 **No setter or deleter, only getter**
210 """
212 if not 'adjacencyMatrix' in self._cache:
213 self._cache['adjacencyMatrix'] = []
214 for i in range(len(self.nl)):
215 self._cache['adjacencyMatrix'].append(self.nl[i].get_connectivity_matrix())
217 return self._cache['adjacencyMatrix']
219 @property
220 def distance_matrix(self):
221 """The distance matrix.
223 If not already done, build a list of distance matrices for all :data:`nl`. See
224 :meth:`ase.neighborlist.get_distance_matrix`.
226 **No setter or deleter, only getter**
227 """
229 if not 'distanceMatrix' in self._cache:
230 self._cache['distanceMatrix'] = []
231 for i in range(len(self.nl)):
232 self._cache['distanceMatrix'].append(get_distance_matrix(self.adjacency_matrix[i]))
234 return self._cache['distanceMatrix']
237 @property
238 def unique_bonds(self):
239 """Get Unique Bonds.
241 :data:`all_bonds` i-j without j-i. This is the upper triangle of the
242 connectivity matrix (i,j), `i < j`
244 """
245 bonds = []
246 for imI in range(len(self.all_bonds)):
247 bonds.append([])
248 for iAtom, bonded in enumerate(self.all_bonds[imI]):
249 bonds[-1].append([ jAtom for jAtom in bonded if jAtom > iAtom ])
251 return bonds
254 def _filter_unique(self, l):
255 """Helper function to filter for unique lists in a list
256 that also contains the reversed items.
257 """
258 r = []
259 #iterate over images
260 for imI in range(len(l)):
261 r.append([])
262 #iterate over atoms
263 for i, tuples in enumerate(l[imI]):
264 #add the ones where i is smaller than the last element
265 r[-1].append([ x for x in tuples if i < x[-1] ])
266 return r
268 def clear_cache(self):
269 """Delete all cached information."""
270 self._cache = {}
272 @property
273 def unique_angles(self):
274 """Get Unique Angles.
276 :data:`all_angles` i-j-k without k-j-i.
278 """
279 return self._filter_unique(self.all_angles)
281 @property
282 def unique_dihedrals(self):
283 """Get Unique Dihedrals.
285 :data:`all_dihedrals` i-j-k-l without l-k-j-i.
287 """
288 return self._filter_unique(self.all_dihedrals)
291 def _get_symbol_idxs(self, imI, sym):
292 """Get list of indices of element *sym*"""
293 if isinstance(imI, int):
294 return [ idx for idx in range(len(self.images[imI])) if self.images[imI][idx].symbol == sym ]
295 else:
296 return [ idx for idx in range(len(imI)) if imI[idx].symbol == sym ]
299 def _idxTuple2SymbolTuple(self, imI, tup):
300 """Converts a tuple of indices to their symbols"""
301 return ( self.images[imI][idx].symbol for idx in tup )
304 def get_bonds(self, A, B, unique=True):
305 """Get bonds from element A to element B.
307 Parameters:
309 A, B: str
310 Get Bonds between elements A and B
311 unique: bool
312 Return the bonds both ways or just one way (A-B and B-A or only A-B)
314 Returns:
316 return: list of lists of tuples
317 return[imageIdx][atomIdx][bondI], each tuple starts with atomIdx.
319 Use :func:`get_values` to convert the returned list to values.
320 """
321 r = []
322 for imI in range(len(self.all_bonds)):
323 r.append([])
324 aIdxs = self._get_symbol_idxs(imI, A)
325 if A != B:
326 bIdxs = self._get_symbol_idxs(imI, B)
327 for idx in aIdxs:
328 bonded = self.all_bonds[imI][idx]
329 if A == B:
330 r[-1].extend([ (idx, x) for x in bonded if ( x in aIdxs ) and ( x > idx ) ])
331 else:
332 r[-1].extend([ (idx, x) for x in bonded if x in bIdxs ])
334 if not unique:
335 r[-1] += [ x[::-1] for x in r[-1] ]
337 return r
340 def get_angles(self, A, B, C, unique=True):
341 """Get angles from given elements A-B-C.
343 Parameters:
345 A, B, C: str
346 Get Angles between elements A, B and C. **B will be the central atom**.
347 unique: bool
348 Return the angles both ways or just one way (A-B-C and C-B-A or only A-B-C)
350 Returns:
352 return: list of lists of tuples
353 return[imageIdx][atomIdx][angleI], each tuple starts with atomIdx.
355 Use :func:`get_values` to convert the returned list to values.
356 """
357 from itertools import product, combinations
358 r = []
359 for imI in range(len(self.all_angles)):
360 r.append([])
361 #Middle Atom is fixed
362 bIdxs = self._get_symbol_idxs(imI, B)
363 for bIdx in bIdxs:
364 bondedA = [ idx for idx in self.all_bonds[imI][bIdx] if self.images[imI][idx].symbol == A ]
365 if len(bondedA) == 0:
366 continue
368 if A != C:
369 bondedC = [ idx for idx in self.all_bonds[imI][bIdx] if self.images[imI][idx].symbol == C ]
370 if len(bondedC) == 0:
371 continue
373 if A == C:
374 extend = [ (x[0], bIdx, x[1]) for x in list(combinations(bondedA, 2)) ]
375 else:
376 extend = list(product(bondedA, [bIdx], bondedC))
378 if not unique:
379 extend += [ x[::-1] for x in extend ]
381 r[-1].extend(extend)
382 return r
385 def get_dihedrals(self, A, B, C, D, unique=True):
386 """Get dihedrals A-B-C-D.
388 Parameters:
390 A, B, C, D: str
391 Get Dihedralss between elements A, B, C and D. **B-C will be the central axis**.
392 unique: bool
393 Return the dihedrals both ways or just one way (A-B-C-D and D-C-B-A or only A-B-C-D)
395 Returns:
397 return: list of lists of tuples
398 return[imageIdx][atomIdx][dihedralI], each tuple starts with atomIdx.
400 Use :func:`get_values` to convert the returned list to values.
401 """
402 r = []
403 for imI in range(len(self.all_dihedrals)):
404 r.append([])
405 #get indices of elements
406 aIdxs = self._get_symbol_idxs(imI, A)
407 bIdxs = self._get_symbol_idxs(imI, B)
408 cIdxs = self._get_symbol_idxs(imI, C)
409 dIdxs = self._get_symbol_idxs(imI, D)
410 for aIdx in aIdxs:
411 dihedrals = [ (aIdx, ) + d for d in self.all_dihedrals[imI][aIdx] if ( d[0] in bIdxs ) and ( d[1] in cIdxs ) and ( d[2] in dIdxs ) ]
412 if not unique:
413 dihedrals += [ d[::-1] for d in dihedrals ]
414 r[-1].extend(dihedrals)
416 return r
419 def get_bond_value(self, imIdx, idxs, mic=True, **kwargs):
420 """Get bond length.
422 Parameters:
424 imIdx: int
425 Index of Image to get value from.
426 idxs: tuple or list of integers
427 Get distance between atoms idxs[0]-idxs[1].
428 mic: bool
429 Passed on to :func:`ase.Atoms.get_distance` for retrieving the value, defaults to True.
430 If the cell of the image is correctly set, there should be no reason to change this.
431 kwargs: options or dict
432 Passed on to :func:`ase.Atoms.get_distance`.
434 Returns:
436 return: float
437 Value returned by image.get_distance.
438 """
439 return self.images[imIdx].get_distance(idxs[0], idxs[1], mic=mic, **kwargs)
441 def get_angle_value(self, imIdx, idxs, mic=True, **kwargs):
442 """Get angle.
444 Parameters:
446 imIdx: int
447 Index of Image to get value from.
448 idxs: tuple or list of integers
449 Get angle between atoms idxs[0]-idxs[1]-idxs[2].
450 mic: bool
451 Passed on to :func:`ase.Atoms.get_angle` for retrieving the value, defaults to True.
452 If the cell of the image is correctly set, there should be no reason to change this.
453 kwargs: options or dict
454 Passed on to :func:`ase.Atoms.get_angle`.
456 Returns:
458 return: float
459 Value returned by image.get_angle.
460 """
461 return self.images[imIdx].get_angle(idxs[0], idxs[1], idxs[2], mic=True, **kwargs)
463 def get_dihedral_value(self, imIdx, idxs, mic=True, **kwargs):
464 """Get dihedral.
466 Parameters:
468 imIdx: int
469 Index of Image to get value from.
470 idxs: tuple or list of integers
471 Get angle between atoms idxs[0]-idxs[1]-idxs[2]-idxs[3].
472 mic: bool
473 Passed on to :func:`ase.Atoms.get_dihedral` for retrieving the value, defaults to True.
474 If the cell of the image is correctly set, there should be no reason to change this.
475 kwargs: options or dict
476 Passed on to :func:`ase.Atoms.get_dihedral`.
478 Returns:
480 return: float
481 Value returned by image.get_dihedral.
482 """
483 return self.images[imIdx].get_dihedral(idxs[0], idxs[1], idxs[2], idxs[3], mic=mic, **kwargs)
485 def get_values(self, inputList, imageIdx=None, mic=True, **kwargs):
486 """Get Bond/Angle/Dihedral values.
488 Parameters:
490 inputList: list of lists of tuples
491 Can be any list provided by :meth:`~ase.geometry.analysis.Analysis.get_bonds`,
492 :meth:`~ase.geometry.analysis.Analysis.get_angles` or
493 :meth:`~ase.geometry.analysis.Analysis.get_dihedrals`.
494 imageIdx: integer or slice
495 The images from :data:`images` to be analyzed. If None, all frames will be analyzed.
496 See :func:`~ase.geometry.analysis.Analysis._get_slice` for details.
497 mic: bool
498 Passed on to :class:`~ase.Atoms` for retrieving the values, defaults to True.
499 If the cells of the images are correctly set, there should be no reason to change this.
500 kwargs: options or dict
501 Passed on to the :class:`~ase.Atoms` classes functions for retrieving the values.
503 Returns:
505 return: list of lists of floats
506 return[imageIdx][valueIdx]. Has the same shape as the *inputList*, instead of each
507 tuple there is a float with the value this tuple yields.
509 The type of value requested is determined from the length of the tuple inputList[0][0].
510 The methods from the :class:`~ase.Atoms` class are used.
511 """
513 sl = self._get_slice(imageIdx)
515 #get method to call from length of inputList
516 if len(inputList[0][0]) == 2:
517 get = self.get_bond_value
518 elif len(inputList[0][0]) == 3:
519 get = self.get_angle_value
520 elif len(inputList[0][0]) == 4:
521 get = self.get_dihedral_value
522 else:
523 raise ValueError("inputList in ase.geometry.analysis.Analysis.get_values has a bad shape.")
525 #check if length of slice and inputList match
526 singleNL = False
527 if len(inputList) != len(self.images[sl]):
528 #only one nl for all images
529 if len(inputList) == 1 and len(self.nl) == 1:
530 singleNL = True
531 else:
532 raise RuntimeError("Length of inputList does not match length of \
533 images requested, but it also is not one item long.")
535 r = []
536 for inputIdx, image in enumerate(self.images[sl]):
537 imageIdx = self.images.index(image)
538 r.append([])
539 #always use first list from input if only a single neighborlist was used
540 if singleNL:
541 inputIdx = 0
542 for tupl in inputList[inputIdx]:
543 r[-1].append(get(imageIdx, tupl, mic=mic, **kwargs))
545 return r
548 def get_rdf(self, rmax, nbins, imageIdx=None, elements=None, return_dists=False):
549 """Get RDF.
551 Wrapper for :meth:`ase.ga.utilities.get_rdf` with more selection possibilities.
553 Parameters:
555 rmax: float
556 Maximum distance of RDF.
557 nbins: int
558 Number of bins to divide RDF.
559 imageIdx: int/slice/None
560 Images to analyze, see :func:`_get_slice` for details.
561 elements: str/int/list/tuple
562 Make partial RDFs.
564 If elements is *None*, a full RDF is calculated. If elements is an *integer* or a *list/tuple
565 of integers*, only those atoms will contribute to the RDF (like a mask). If elements
566 is a *string* or a *list/tuple of strings*, only Atoms of those elements will contribute.
568 Returns:
570 return: list of lists / list of tuples of lists
571 If return_dists is True, the returned tuples contain (rdf, distances). Otherwise
572 only rdfs for each image are returned.
573 """
575 sl = self._get_slice(imageIdx)
577 r = []
578 el = None
580 for image in self.images[sl]:
581 if elements is None:
582 tmpImage = image
583 #integers
584 elif isinstance(elements, int):
585 tmpImage = Atoms(cell=image.get_cell(), pbc=image.get_pbc())
586 tmpImage.append(image[elements])
587 #strings
588 elif isinstance(elements, str):
589 tmpImage = Atoms(cell=image.get_cell(), pbc=image.get_pbc())
590 for idx in self._get_symbol_idxs(image, elements):
591 tmpImage.append(image[idx])
592 #lists
593 elif isinstance(elements, list) or isinstance(elements, tuple):
594 #list of ints
595 if all(isinstance(x, int) for x in elements):
596 if len(elements) == 2:
597 #use builtin get_rdf mask
598 el = elements
599 tmpImage = image
600 else:
601 #create dummy image
602 tmpImage = Atoms(cell=image.get_cell(), pbc=image.get_pbc())
603 for idx in elements:
604 tmpImage.append(image[idx])
605 #list of strings
606 elif all(isinstance(x, str) for x in elements):
607 tmpImage = Atoms(cell=image.get_cell(), pbc=image.get_pbc())
608 for element in elements:
609 for idx in self._get_symbol_idxs(image, element):
610 tmpImage.append(image[idx])
611 else:
612 raise ValueError("Unsupported type of elements given in ase.geometry.analysis.Analysis.get_rdf!")
613 else:
614 raise ValueError("Unsupported type of elements given in ase.geometry.analysis.Analysis.get_rdf!")
616 r.append(get_rdf(tmpImage, rmax, nbins, elements=el, no_dists=(not return_dists)))
617 return r