Coverage for ase / geometry / analysis.py: 64.29%
266 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3# flake8: noqa
4"""Tools for analyzing instances of :class:`~ase.Atoms`
5"""
7from typing import List
9import numpy as np
11from ase import Atoms
12from ase.geometry.rdf import get_containing_cell_length, get_rdf
13from ase.neighborlist import (build_neighbor_list, get_distance_indices,
14 get_distance_matrix)
15from ase.utils import deprecated
17__all__ = ['Analysis']
20def get_max_containing_cell_length(images: List[Atoms]):
21 i2diff = np.zeros(3)
22 for image in images:
23 np.maximum(get_containing_cell_length(image), i2diff, out=i2diff)
24 return i2diff
27def get_max_volume_estimate(images: List[Atoms]) -> float:
28 return np.prod(get_max_containing_cell_length(images))
31class Analysis:
32 """Analysis class
34 Parameters for initialization:
36 images: :class:`~ase.Atoms` object or list of such
37 Images to analyze.
38 nl: None, :class:`~ase.neighborlist.NeighborList` object or list of such
39 Neighborlist(s) for the given images. One or nImages, depending if bonding
40 pattern changes or is constant. Using one Neigborlist greatly improves speed.
41 kwargs: options, dict
42 Arguments for constructing :class:`~ase.neighborlist.NeighborList` object if :data:`nl` is None.
44 The choice of ``bothways=True`` for the :class:`~ase.neighborlist.NeighborList` object
45 will not influence the amount of bonds/angles/dihedrals you get, all are reported
46 in both directions. Use the *unique*-labeled properties to get lists without
47 duplicates.
48 """
50 def __init__(self, images, nl=None, **kwargs):
51 self.images = images
53 if isinstance(nl, list):
54 assert len(nl) == self.nImages
55 self._nl = nl
56 elif nl is not None:
57 self._nl = [nl]
58 else:
59 self._nl = [build_neighbor_list(self.images[0], **kwargs)]
61 self._cache = {}
63 def _get_slice(self, imageIdx):
64 """Return a slice from user input.
65 Using *imageIdx* (can be integer or slice) the analyzed frames can be specified.
66 If *imageIdx* is None, all frames will be analyzed.
67 """
68 # get slice from imageIdx
69 if isinstance(imageIdx, int):
70 sl = slice(imageIdx, imageIdx + 1)
71 elif isinstance(imageIdx, slice):
72 sl = imageIdx
73 elif imageIdx is None:
74 sl = slice(0, None)
75 else:
76 raise ValueError(
77 "Unsupported type for imageIdx in ase.geometry.analysis.Analysis._get_slice")
78 return sl
80 @property
81 def images(self):
82 """Images.
84 Set during initialization but can also be set later.
85 """
86 return self._images
88 @images.setter
89 def images(self, images):
90 """Set images"""
91 if isinstance(images, list):
92 self._images = images
93 else:
94 self._images = [images]
96 @images.deleter
97 def images(self):
98 """Delete images"""
99 self._images = None
101 @property
102 def nImages(self):
103 """Number of Images in this instance.
105 Cannot be set, is determined automatically.
106 """
107 return len(self.images)
109 @property
110 def nl(self):
111 """Neighbor Lists in this instance.
113 Set during initialization.
115 **No setter or deleter, only getter**
116 """
117 return self._nl
119 def _get_all_x(self, distance):
120 """Helper function to get bonds, angles, dihedrals"""
121 maxIter = self.nImages
122 if len(self.nl) == 1:
123 maxIter = 1
125 xList = []
126 for i in range(maxIter):
127 xList.append(get_distance_indices(
128 self.distance_matrix[i], distance))
130 return xList
132 @property
133 def all_bonds(self):
134 """All Bonds.
136 A list with indices of bonded atoms for each neighborlist in *self*.
137 Atom i is connected to all atoms inside result[i]. Duplicates from PBCs are
138 removed. See also :data:`unique_bonds`.
140 **No setter or deleter, only getter**
141 """
142 if 'allBonds' not in self._cache:
143 self._cache['allBonds'] = self._get_all_x(1)
145 return self._cache['allBonds']
147 @property
148 def all_angles(self):
149 """All angles
151 A list with indices of atoms in angles for each neighborlist in *self*.
152 Atom i forms an angle to the atoms inside the tuples in result[i]:
153 i -- result[i][x][0] -- result[i][x][1]
154 where x is in range(number of angles from i). See also :data:`unique_angles`.
156 **No setter or deleter, only getter**
157 """
158 if 'allAngles' not in self._cache:
159 self._cache['allAngles'] = []
160 distList = self._get_all_x(2)
162 for imI in range(len(distList)):
163 self._cache['allAngles'].append([])
164 # iterate over second neighbors of all atoms
165 for iAtom, secNeighs in enumerate(distList[imI]):
166 self._cache['allAngles'][-1].append([])
167 if len(secNeighs) == 0:
168 continue
169 firstNeighs = self.all_bonds[imI][iAtom]
170 # iterate over second neighbors of iAtom
171 for kAtom in secNeighs:
172 relevantFirstNeighs = [
173 idx for idx in firstNeighs if kAtom in self.all_bonds[imI][idx]]
174 # iterate over all atoms that are connected to iAtom and kAtom
175 for jAtom in relevantFirstNeighs:
176 self._cache['allAngles'][-1][-1].append(
177 (jAtom, kAtom))
179 return self._cache['allAngles']
181 @property
182 def all_dihedrals(self):
183 """All dihedrals
185 Returns a list with indices of atoms in dihedrals for each neighborlist in this instance.
186 Atom i forms a dihedral to the atoms inside the tuples in result[i]:
187 i -- result[i][x][0] -- result[i][x][1] -- result[i][x][2]
188 where x is in range(number of dihedrals from i). See also :data:`unique_dihedrals`.
190 **No setter or deleter, only getter**
191 """
192 if 'allDihedrals' not in self._cache:
193 self._cache['allDihedrals'] = []
194 distList = self._get_all_x(3)
196 for imI in range(len(distList)):
197 self._cache['allDihedrals'].append([])
198 for iAtom, thirdNeighs in enumerate(distList[imI]):
199 self._cache['allDihedrals'][-1].append([])
200 if len(thirdNeighs) == 0:
201 continue
202 anglesI = self.all_angles[imI][iAtom]
203 # iterate over third neighbors of iAtom
204 for lAtom in thirdNeighs:
205 secondNeighs = [angle[-1] for angle in anglesI]
206 firstNeighs = [angle[0] for angle in anglesI]
207 relevantSecondNeighs = [
208 idx for idx in secondNeighs if lAtom in self.all_bonds[imI][idx]]
209 relevantFirstNeighs = [
210 firstNeighs[secondNeighs.index(idx)] for idx in relevantSecondNeighs]
211 # iterate over all atoms that are connected to iAtom and lAtom
212 for jAtom, kAtom in zip(relevantFirstNeighs, relevantSecondNeighs):
213 # remove dihedrals in circles
214 tupl = (jAtom, kAtom, lAtom)
215 if len(set((iAtom, ) + tupl)) != 4:
216 continue
217 # avoid duplicates
218 elif tupl in self._cache['allDihedrals'][-1][-1]:
219 continue
220 elif iAtom in tupl:
221 raise RuntimeError(
222 "Something is wrong in analysis.all_dihedrals!")
223 self._cache['allDihedrals'][-1][-1].append(
224 (jAtom, kAtom, lAtom))
226 return self._cache['allDihedrals']
228 @property
229 def adjacency_matrix(self):
230 """The adjacency/connectivity matrix.
232 If not already done, build a list of adjacency matrices for all :data:`nl`.
234 **No setter or deleter, only getter**
235 """
237 if 'adjacencyMatrix' not in self._cache:
238 self._cache['adjacencyMatrix'] = []
239 for i in range(len(self.nl)):
240 self._cache['adjacencyMatrix'].append(
241 self.nl[i].get_connectivity_matrix())
243 return self._cache['adjacencyMatrix']
245 @property
246 def distance_matrix(self):
247 """The distance matrix.
249 If not already done, build a list of distance matrices for all :data:`nl`. See
250 :meth:`ase.neighborlist.get_distance_matrix`.
252 **No setter or deleter, only getter**
253 """
255 if 'distanceMatrix' not in self._cache:
256 self._cache['distanceMatrix'] = []
257 for i in range(len(self.nl)):
258 self._cache['distanceMatrix'].append(
259 get_distance_matrix(self.adjacency_matrix[i]))
261 return self._cache['distanceMatrix']
263 @property
264 def unique_bonds(self):
265 """Get Unique Bonds.
267 :data:`all_bonds` i-j without j-i. This is the upper triangle of the
268 connectivity matrix (i,j), `i < j`
270 """
271 bonds = []
272 for imI in range(len(self.all_bonds)):
273 bonds.append([])
274 for iAtom, bonded in enumerate(self.all_bonds[imI]):
275 bonds[-1].append([jAtom for jAtom in bonded if jAtom > iAtom])
277 return bonds
279 def _filter_unique(self, l):
280 """Helper function to filter for unique lists in a list
281 that also contains the reversed items.
282 """
283 r = []
284 # iterate over images
285 for imI in range(len(l)):
286 r.append([])
287 # iterate over atoms
288 for i, tuples in enumerate(l[imI]):
289 # add the ones where i is smaller than the last element
290 r[-1].append([x for x in tuples if i < x[-1]])
291 return r
293 def clear_cache(self):
294 """Delete all cached information."""
295 self._cache = {}
297 @property
298 def unique_angles(self):
299 """Get Unique Angles.
301 :data:`all_angles` i-j-k without k-j-i.
303 """
304 return self._filter_unique(self.all_angles)
306 @property
307 def unique_dihedrals(self):
308 """Get Unique Dihedrals.
310 :data:`all_dihedrals` i-j-k-l without l-k-j-i.
312 """
313 return self._filter_unique(self.all_dihedrals)
315 def _get_symbol_idxs(self, imI, sym):
316 """Get list of indices of element *sym*"""
317 if isinstance(imI, int):
318 return [idx for idx in range(len(self.images[imI])) if self.images[imI][idx].symbol == sym]
319 else:
320 return [idx for idx in range(len(imI)) if imI[idx].symbol == sym]
322 def _idxTuple2SymbolTuple(self, imI, tup):
323 """Converts a tuple of indices to their symbols"""
324 return (self.images[imI][idx].symbol for idx in tup)
326 def get_bonds(self, A, B, unique=True):
327 """Get bonds from element A to element B.
329 Parameters:
331 A, B: str
332 Get Bonds between elements A and B
333 unique: bool
334 Return the bonds both ways or just one way (A-B and B-A or only A-B)
336 Returns:
338 return: list of lists of tuples
339 return[imageIdx][atomIdx][bondI], each tuple starts with atomIdx.
341 Use :func:`get_values` to convert the returned list to values.
342 """
343 r = []
344 for imI in range(len(self.all_bonds)):
345 r.append([])
346 aIdxs = self._get_symbol_idxs(imI, A)
347 if A != B:
348 bIdxs = self._get_symbol_idxs(imI, B)
349 for idx in aIdxs:
350 bonded = self.all_bonds[imI][idx]
351 if A == B:
352 r[-1].extend([(idx, x)
353 for x in bonded if (x in aIdxs) and (x > idx)])
354 else:
355 r[-1].extend([(idx, x) for x in bonded if x in bIdxs])
357 if not unique:
358 r[-1] += [x[::-1] for x in r[-1]]
360 return r
362 def get_angles(self, A, B, C, unique=True):
363 """Get angles from given elements A-B-C.
365 Parameters:
367 A, B, C: str
368 Get Angles between elements A, B and C. **B will be the central atom**.
369 unique: bool
370 Return the angles both ways or just one way (A-B-C and C-B-A or only A-B-C)
372 Returns:
374 return: list of lists of tuples
375 return[imageIdx][atomIdx][angleI], each tuple starts with atomIdx.
377 Use :func:`get_values` to convert the returned list to values.
378 """
379 from itertools import combinations, product
380 r = []
381 for imI in range(len(self.all_angles)):
382 r.append([])
383 # Middle Atom is fixed
384 bIdxs = self._get_symbol_idxs(imI, B)
385 for bIdx in bIdxs:
386 bondedA = [idx for idx in self.all_bonds[imI]
387 [bIdx] if self.images[imI][idx].symbol == A]
388 if len(bondedA) == 0:
389 continue
391 if A != C:
392 bondedC = [idx for idx in self.all_bonds[imI]
393 [bIdx] if self.images[imI][idx].symbol == C]
394 if len(bondedC) == 0:
395 continue
397 if A == C:
398 extend = [(x[0], bIdx, x[1])
399 for x in list(combinations(bondedA, 2))]
400 else:
401 extend = list(product(bondedA, [bIdx], bondedC))
403 if not unique:
404 extend += [x[::-1] for x in extend]
406 r[-1].extend(extend)
407 return r
409 def get_dihedrals(self, A, B, C, D, unique=True):
410 """Get dihedrals A-B-C-D.
412 Parameters:
414 A, B, C, D: str
415 Get Dihedralss between elements A, B, C and D. **B-C will be the central axis**.
416 unique: bool
417 Return the dihedrals both ways or just one way (A-B-C-D and D-C-B-A or only A-B-C-D)
419 Returns:
421 return: list of lists of tuples
422 return[imageIdx][atomIdx][dihedralI], each tuple starts with atomIdx.
424 Use :func:`get_values` to convert the returned list to values.
425 """
426 r = []
427 for imI in range(len(self.all_dihedrals)):
428 r.append([])
429 # get indices of elements
430 aIdxs = self._get_symbol_idxs(imI, A)
431 bIdxs = self._get_symbol_idxs(imI, B)
432 cIdxs = self._get_symbol_idxs(imI, C)
433 dIdxs = self._get_symbol_idxs(imI, D)
434 for aIdx in aIdxs:
435 dihedrals = [(aIdx, ) + d for d in self.all_dihedrals[imI][aIdx]
436 if (d[0] in bIdxs) and (d[1] in cIdxs) and (d[2] in dIdxs)]
437 if not unique:
438 dihedrals += [d[::-1] for d in dihedrals]
439 r[-1].extend(dihedrals)
441 return r
443 def get_bond_value(self, imIdx, idxs, mic=True, **kwargs):
444 """Get bond length.
446 Parameters:
448 imIdx: int
449 Index of Image to get value from.
450 idxs: tuple or list of integers
451 Get distance between atoms idxs[0]-idxs[1].
452 mic: bool
453 Passed on to :func:`ase.Atoms.get_distance` for retrieving the value, defaults to True.
454 If the cell of the image is correctly set, there should be no reason to change this.
455 kwargs: options or dict
456 Passed on to :func:`ase.Atoms.get_distance`.
458 Returns:
460 return: float
461 Value returned by image.get_distance.
462 """
463 return self.images[imIdx].get_distance(idxs[0], idxs[1], mic=mic, **kwargs)
465 def get_angle_value(self, imIdx, idxs, mic=True, **kwargs):
466 """Get angle.
468 Parameters:
470 imIdx: int
471 Index of Image to get value from.
472 idxs: tuple or list of integers
473 Get angle between atoms idxs[0]-idxs[1]-idxs[2].
474 mic: bool
475 Passed on to :func:`ase.Atoms.get_angle` for retrieving the value, defaults to True.
476 If the cell of the image is correctly set, there should be no reason to change this.
477 kwargs: options or dict
478 Passed on to :func:`ase.Atoms.get_angle`.
480 Returns:
482 return: float
483 Value returned by image.get_angle.
484 """
485 return self.images[imIdx].get_angle(idxs[0], idxs[1], idxs[2], mic=True, **kwargs)
487 def get_dihedral_value(self, imIdx, idxs, mic=True, **kwargs):
488 """Get dihedral.
490 Parameters:
492 imIdx: int
493 Index of Image to get value from.
494 idxs: tuple or list of integers
495 Get angle between atoms idxs[0]-idxs[1]-idxs[2]-idxs[3].
496 mic: bool
497 Passed on to :func:`ase.Atoms.get_dihedral` for retrieving the value, defaults to True.
498 If the cell of the image is correctly set, there should be no reason to change this.
499 kwargs: options or dict
500 Passed on to :func:`ase.Atoms.get_dihedral`.
502 Returns:
504 return: float
505 Value returned by image.get_dihedral.
506 """
507 return self.images[imIdx].get_dihedral(idxs[0], idxs[1], idxs[2], idxs[3], mic=mic, **kwargs)
509 def get_values(self, inputList, imageIdx=None, mic=True, **kwargs):
510 """Get Bond/Angle/Dihedral values.
512 Parameters:
514 inputList: list of lists of tuples
515 Can be any list provided by :meth:`~ase.geometry.analysis.Analysis.get_bonds`,
516 :meth:`~ase.geometry.analysis.Analysis.get_angles` or
517 :meth:`~ase.geometry.analysis.Analysis.get_dihedrals`.
518 imageIdx: integer or slice
519 The images from :data:`images` to be analyzed. If None, all frames will be analyzed.
520 See :func:`~ase.geometry.analysis.Analysis._get_slice` for details.
521 mic: bool
522 Passed on to :class:`~ase.Atoms` for retrieving the values, defaults to True.
523 If the cells of the images are correctly set, there should be no reason to change this.
524 kwargs: options or dict
525 Passed on to the :class:`~ase.Atoms` classes functions for retrieving the values.
527 Returns:
529 return: list of lists of floats
530 return[imageIdx][valueIdx]. Has the same shape as the *inputList*, instead of each
531 tuple there is a float with the value this tuple yields.
533 The type of value requested is determined from the length of the tuple inputList[0][0].
534 The methods from the :class:`~ase.Atoms` class are used.
535 """
537 sl = self._get_slice(imageIdx)
539 # get method to call from length of inputList
540 if len(inputList[0][0]) == 2:
541 get = self.get_bond_value
542 elif len(inputList[0][0]) == 3:
543 get = self.get_angle_value
544 elif len(inputList[0][0]) == 4:
545 get = self.get_dihedral_value
546 else:
547 raise ValueError(
548 "inputList in ase.geometry.analysis.Analysis.get_values has a bad shape.")
550 # check if length of slice and inputList match
551 singleNL = False
552 if len(inputList) != len(self.images[sl]):
553 # only one nl for all images
554 if len(inputList) == 1 and len(self.nl) == 1:
555 singleNL = True
556 else:
557 raise RuntimeError("Length of inputList does not match length of \
558 images requested, but it also is not one item long.")
560 r = []
561 for inputIdx, image in enumerate(self.images[sl]):
562 imageIdx = self.images.index(image)
563 r.append([])
564 # always use first list from input if only a single neighborlist was used
565 if singleNL:
566 inputIdx = 0
567 for tupl in inputList[inputIdx]:
568 r[-1].append(get(imageIdx, tupl, mic=mic, **kwargs))
570 return r
572 def get_max_volume_estimate(self):
573 return get_max_volume_estimate(self.images)
575 @deprecated("Use `ase.geometry.rdf.get_rdf` instead.")
576 def get_rdf(self, rmax, nbins, imageIdx=None, elements=None, return_dists=False,
577 volume: float | None = None):
578 """Get RDF.
580 Wrapper for :func:`ase.geometry.rdf.get_rdf` with more selection possibilities.
582 .. deprecated:: 3.28.0
583 This method has bugs when ``elements`` is not ``None`` and will be
584 removed soon. Use :func:`ase.geometry.rdf.get_rdf` instead.
586 Parameters:
588 rmax: float
589 Maximum distance of RDF.
590 nbins: int
591 Number of bins to divide RDF.
592 imageIdx: int/slice/None
593 Images to analyze, see :func:`_get_slice` for details.
594 elements: str/int/list/tuple
595 Make partial RDFs.
597 If elements is *None*, a full RDF is calculated. If elements is an *integer* or a *list/tuple
598 of integers*, only those atoms will contribute to the RDF (like a mask). If elements
599 is a *string* or a *list/tuple of strings*, only Atoms of those elements will contribute.
601 Returns:
603 return: list of lists / list of tuples of lists
604 If return_dists is True, the returned tuples contain (rdf, distances). Otherwise
605 only rdfs for each image are returned.
606 """
608 if elements is not None:
609 raise RuntimeError(
610 'There are bugs for `elements` is not `None`. '
611 'Use `ase.geometry.rdf.get_rdf` instead.'
612 )
614 sl = self._get_slice(imageIdx)
616 ls_rdf = []
617 el = None
619 for image in self.images[sl]:
620 if elements is None:
621 tmp_image = image
622 # integers
623 elif isinstance(elements, int):
624 tmp_image = Atoms(cell=image.get_cell(), pbc=image.get_pbc())
625 tmp_image.append(image[elements])
626 # strings
627 elif isinstance(elements, str):
628 tmp_image = Atoms(cell=image.get_cell(), pbc=image.get_pbc())
629 for idx in self._get_symbol_idxs(image, elements):
630 tmp_image.append(image[idx])
631 # lists
632 elif isinstance(elements, (list, tuple)):
633 # list of ints
634 if all(isinstance(x, int) for x in elements):
635 if len(elements) == 2:
636 # use builtin get_rdf mask
637 el = elements
638 tmp_image = image
639 else:
640 # create dummy image
641 tmp_image = Atoms(
642 cell=image.get_cell(), pbc=image.get_pbc())
643 for idx in elements:
644 tmp_image.append(image[idx])
645 # list of strings
646 elif all(isinstance(x, str) for x in elements):
647 tmp_image = Atoms(cell=image.get_cell(),
648 pbc=image.get_pbc())
649 for element in elements:
650 for idx in self._get_symbol_idxs(image, element):
651 tmp_image.append(image[idx])
652 else:
653 raise ValueError(
654 "Unsupported type of elements given in ase.geometry.analysis.Analysis.get_rdf!")
655 else:
656 raise ValueError(
657 "Unsupported type of elements given in ase.geometry.analysis.Analysis.get_rdf!")
659 rdf = get_rdf(tmp_image, rmax, nbins, elements=el, no_dists=(not return_dists),
660 volume=volume)
661 ls_rdf.append(rdf)
663 return ls_rdf