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

1import math 

2from typing import List, Optional, Tuple, Union 

3 

4import numpy as np 

5 

6from ase import Atoms 

7from ase.cell import Cell 

8 

9 

10class CellTooSmall(Exception): 

11 pass 

12 

13 

14class VolumeNotDefined(Exception): 

15 pass 

16 

17 

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. 

26 

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. 

31 

32 Parameters: 

33 

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. 

38 

39 nbins : int 

40 Number of bins to divide the rdf into. 

41 

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. 

46 

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. 

50 

51 no_dists : bool 

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

53 

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 """ 

58 

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 

63 

64 check_cell_and_r_max(atoms, rmax) 

65 

66 dm = distance_matrix 

67 if dm is None: 

68 dm = atoms.get_all_distances(mic=True) 

69 

70 rdf = np.zeros(nbins + 1) 

71 dr = float(rmax / nbins) 

72 

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

74 natoms = len(atoms) 

75 

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) 

80 

81 indices_triu = np.triu(indices) 

82 for index in range(nbins + 1): 

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

84 

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 

89 

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 

95 

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

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

98 

99 if no_dists: 

100 return rdf[1:] 

101 

102 return rdf[1:], rr 

103 

104 

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

106 cell = atoms.get_cell() 

107 pbc = atoms.get_pbc() 

108 

109 vol = atoms.cell.volume 

110 

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}') 

121 

122 

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 

134 

135 

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 

139 

140 

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

142 return np.prod(get_containing_cell_length(atoms))