Coverage for /builds/kinetik161/ase/ase/geometry/rdf.py: 100.00%
66 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import math
2from typing import List, Optional, Tuple, Union
4import numpy as np
6from ase import Atoms
7from ase.cell import Cell
10class CellTooSmall(Exception):
11 pass
14class VolumeNotDefined(Exception):
15 pass
18def get_rdf(atoms: Atoms, rmax: float, nbins: int,
19 distance_matrix: Optional[np.ndarray] = None,
20 elements: Optional[Union[List[int], Tuple]] = None,
21 no_dists: Optional[bool] = False,
22 volume: Optional[float] = None):
23 """Returns two numpy arrays; the radial distribution function
24 and the corresponding distances of the supplied atoms object.
25 If no_dists = True then only the first array is returned.
27 Note that the rdf is computed following the standard solid state
28 definition which uses the cell volume in the normalization.
29 This may or may not be appropriate in cases where one or more
30 directions is non-periodic.
32 Parameters:
34 rmax : float
35 The maximum distance that will contribute to the rdf.
36 The unit cell should be large enough so that it encloses a
37 sphere with radius rmax in the periodic directions.
39 nbins : int
40 Number of bins to divide the rdf into.
42 distance_matrix : numpy.array
43 An array of distances between atoms, typically
44 obtained by atoms.get_all_distances().
45 Default None meaning that it will be calculated.
47 elements : list or tuple
48 List of two atomic numbers. If elements is not None the partial
49 rdf for the supplied elements will be returned.
51 no_dists : bool
52 If True then the second array with rdf distances will not be returned.
54 volume : float or None
55 Optionally specify the volume of the system. If specified, the volume
56 will be used instead atoms.cell.volume.
57 """
59 # First check whether the cell is sufficiently large
60 vol = atoms.cell.volume if volume is None else volume
61 if vol < 1.0e-10:
62 raise VolumeNotDefined
64 check_cell_and_r_max(atoms, rmax)
66 dm = distance_matrix
67 if dm is None:
68 dm = atoms.get_all_distances(mic=True)
70 rdf = np.zeros(nbins + 1)
71 dr = float(rmax / nbins)
73 indices = np.asarray(np.ceil(dm / dr), dtype=int)
74 natoms = len(atoms)
76 if elements is None:
77 # Coefficients to use for normalization
78 phi = natoms / vol
79 norm = 2.0 * math.pi * dr * phi * len(atoms)
81 indices_triu = np.triu(indices)
82 for index in range(nbins + 1):
83 rdf[index] = np.count_nonzero(indices_triu == index)
85 else:
86 i_indices = np.where(atoms.numbers == elements[0])[0]
87 phi = len(i_indices) / vol
88 norm = 4.0 * math.pi * dr * phi * natoms
90 for i in i_indices:
91 for j in np.where(atoms.numbers == elements[1])[0]:
92 index = indices[i, j]
93 if index <= nbins:
94 rdf[index] += 1
96 rr = np.arange(dr / 2, rmax, dr)
97 rdf[1:] /= norm * (rr * rr + (dr * dr / 12))
99 if no_dists:
100 return rdf[1:]
102 return rdf[1:], rr
105def check_cell_and_r_max(atoms: Atoms, rmax: float) -> None:
106 cell = atoms.get_cell()
107 pbc = atoms.get_pbc()
109 vol = atoms.cell.volume
111 for i in range(3):
112 if pbc[i]:
113 axb = np.cross(cell[(i + 1) % 3, :], cell[(i + 2) % 3, :])
114 h = vol / np.linalg.norm(axb)
115 if h < 2 * rmax:
116 recommended_r_max = get_recommended_r_max(cell, pbc)
117 raise CellTooSmall(
118 'The cell is not large enough in '
119 f'direction {i}: {h:.3f} < 2*rmax={2 * rmax: .3f}. '
120 f'Recommended rmax = {recommended_r_max}')
123def get_recommended_r_max(cell: Cell, pbc: List[bool]) -> float:
124 recommended_r_max = 5.0
125 vol = cell.volume
126 for i in range(3):
127 if pbc[i]:
128 axb = np.cross(cell[(i + 1) % 3, :], # type: ignore[index]
129 cell[(i + 2) % 3, :]) # type: ignore[index]
130 h = vol / np.linalg.norm(axb)
131 assert isinstance(h, float) # mypy
132 recommended_r_max = min(h / 2 * 0.99, recommended_r_max)
133 return recommended_r_max
136def get_containing_cell_length(atoms: Atoms) -> np.ndarray:
137 atom2xyz = atoms.get_positions()
138 return np.amax(atom2xyz, axis=0) - np.amin(atom2xyz, axis=0) + 2.0
141def get_volume_estimate(atoms: Atoms) -> float:
142 return np.prod(get_containing_cell_length(atoms))