Coverage for /builds/ase/ase/ase/geometry/rdf.py: 100.00%

1# fmt: off 

2 

3import math 

4from typing import List, Optional, Tuple, Union 

5 

6import numpy as np 

7 

8from ase import Atoms 

9from ase.cell import Cell 

10 

11 

12class CellTooSmall(Exception): 

13 pass 

14 

15 

16class VolumeNotDefined(Exception): 

17 pass 

18 

19 

20def get_rdf(atoms: Atoms, rmax: float, nbins: int, 

21 distance_matrix: Optional[np.ndarray] = None, 

22 elements: Optional[Union[List[int], Tuple]] = None, 

23 no_dists: Optional[bool] = False, 

24 volume: Optional[float] = None): 

25 """Returns two numpy arrays; the radial distribution function 

26 and the corresponding distances of the supplied atoms object. 

27 If no_dists = True then only the first array is returned. 

28 

29 Note that the rdf is computed following the standard solid state 

30 definition which uses the cell volume in the normalization. 

31 This may or may not be appropriate in cases where one or more 

32 directions is non-periodic. 

33 

34 Parameters: 

35 

36 rmax : float 

37 The maximum distance that will contribute to the rdf. 

38 The unit cell should be large enough so that it encloses a 

39 sphere with radius rmax in the periodic directions. 

40 

41 nbins : int 

42 Number of bins to divide the rdf into. 

43 

44 distance_matrix : numpy.array 

45 An array of distances between atoms, typically 

46 obtained by atoms.get_all_distances(). 

47 Default None meaning that it will be calculated. 

48 

49 elements : list or tuple 

50 List of two atomic numbers. If elements is not None the partial 

51 rdf for the supplied elements will be returned. 

52 

53 no_dists : bool 

54 If True then the second array with rdf distances will not be returned. 

55 

56 volume : float or None 

57 Optionally specify the volume of the system. If specified, the volume 

58 will be used instead atoms.cell.volume. 

59 """ 

60 

61 # First check whether the cell is sufficiently large 

62 vol = atoms.cell.volume if volume is None else volume 

63 if vol < 1.0e-10: 

64 raise VolumeNotDefined 

65 

66 check_cell_and_r_max(atoms, rmax) 

67 

68 dm = distance_matrix 

69 if dm is None: 

70 dm = atoms.get_all_distances(mic=True) 

71 

72 rdf = np.zeros(nbins + 1) 

73 dr = float(rmax / nbins) 

74 

75 indices = np.asarray(np.ceil(dm / dr), dtype=int) 

76 natoms = len(atoms) 

77 

78 if elements is None: 

79 # Coefficients to use for normalization 

80 phi = natoms / vol 

81 norm = 2.0 * math.pi * dr * phi * len(atoms) 

82 

83 indices_triu = np.triu(indices) 

84 for index in range(nbins + 1): 

85 rdf[index] = np.count_nonzero(indices_triu == index) 

86 

87 else: 

88 i_indices = np.where(atoms.numbers == elements[0])[0] 

89 phi = len(i_indices) / vol 

90 norm = 4.0 * math.pi * dr * phi * natoms 

91 

92 for i in i_indices: 

93 for j in np.where(atoms.numbers == elements[1])[0]: 

94 index = indices[i, j] 

95 if index <= nbins: 

96 rdf[index] += 1 

97 

98 rr = np.arange(dr / 2, rmax, dr) 

99 rdf[1:] /= norm * (rr * rr + (dr * dr / 12)) 

100 

101 if no_dists: 

102 return rdf[1:] 

103 

104 return rdf[1:], rr 

105 

106 

107def check_cell_and_r_max(atoms: Atoms, rmax: float) -> None: 

108 cell = atoms.get_cell() 

109 pbc = atoms.get_pbc() 

110 

111 vol = atoms.cell.volume 

112 

113 for i in range(3): 

114 if pbc[i]: 

115 axb = np.cross(cell[(i + 1) % 3, :], cell[(i + 2) % 3, :]) 

116 h = vol / np.linalg.norm(axb) 

117 if h < 2 * rmax: 

118 recommended_r_max = get_recommended_r_max(cell, pbc) 

119 raise CellTooSmall( 

120 'The cell is not large enough in ' 

121 f'direction {i}: {h:.3f} < 2*rmax={2 * rmax: .3f}. ' 

122 f'Recommended rmax = {recommended_r_max}') 

123 

124 

125def get_recommended_r_max(cell: Cell, pbc: List[bool]) -> float: 

126 recommended_r_max = 5.0 

127 vol = cell.volume 

128 for i in range(3): 

129 if pbc[i]: 

130 axb = np.cross(cell[(i + 1) % 3, :], # type: ignore[index] 

131 cell[(i + 2) % 3, :]) # type: ignore[index] 

132 h = vol / np.linalg.norm(axb) 

133 assert isinstance(h, float) # mypy 

134 recommended_r_max = min(h / 2 * 0.99, recommended_r_max) 

135 return recommended_r_max 

136 

137 

138def get_containing_cell_length(atoms: Atoms) -> np.ndarray: 

139 atom2xyz = atoms.get_positions() 

140 return np.amax(atom2xyz, axis=0) - np.amin(atom2xyz, axis=0) + 2.0 

141 

142 

143def get_volume_estimate(atoms: Atoms) -> float: 

144 return np.prod(get_containing_cell_length(atoms))