1import math 

2from typing import List, Optional, Tuple, Union 

3 

4import numpy as np 

5from ase import Atoms 

6from ase.cell import Cell 

7 

8 

9class CellTooSmall(Exception): 

10 pass 

11 

12 

13class VolumeNotDefined(Exception): 

14 pass 

15 

16 

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

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

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

20 no_dists: Optional[bool] = False, 

21 volume: Optional[float] = None): 

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

23 and the corresponding distances of the supplied atoms object. 

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

25 

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

27 definition which uses the cell volume in the normalization. 

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

29 directions is non-periodic. 

30 

31 Parameters: 

32 

33 rmax : float 

34 The maximum distance that will contribute to the rdf. 

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

36 sphere with radius rmax in the periodic directions. 

37 

38 nbins : int 

39 Number of bins to divide the rdf into. 

40 

41 distance_matrix : numpy.array 

42 An array of distances between atoms, typically 

43 obtained by atoms.get_all_distances(). 

44 Default None meaning that it will be calculated. 

45 

46 elements : list or tuple 

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

48 rdf for the supplied elements will be returned. 

49 

50 no_dists : bool 

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

52 

53 volume : float or None 

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

55 will be used instead atoms.cell.volume. 

56 """ 

57 

58 # First check whether the cell is sufficiently large 

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

60 if vol < 1.0e-10: 

61 raise VolumeNotDefined 

62 

63 check_cell_and_r_max(atoms, rmax) 

64 

65 dm = distance_matrix 

66 if dm is None: 

67 dm = atoms.get_all_distances(mic=True) 

68 

69 rdf = np.zeros(nbins + 1) 

70 dr = float(rmax / nbins) 

71 

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

73 natoms = len(atoms) 

74 

75 if elements is None: 

76 # Coefficients to use for normalization 

77 phi = natoms / vol 

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

79 

80 indices_triu = np.triu(indices) 

81 for index in range(nbins + 1): 

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

83 

84 else: 

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

86 phi = len(i_indices) / vol 

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

88 

89 for i in i_indices: 

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

91 index = indices[i, j] 

92 if index <= nbins: 

93 rdf[index] += 1 

94 

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

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

97 

98 if no_dists: 

99 return rdf[1:] 

100 

101 return rdf[1:], rr 

102 

103 

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

105 cell = atoms.get_cell() 

106 pbc = atoms.get_pbc() 

107 

108 vol = atoms.cell.volume 

109 

110 for i in range(3): 

111 if pbc[i]: 

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

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

114 if h < 2 * rmax: 

115 recommended_r_max = get_recommended_r_max(cell, pbc) 

116 raise CellTooSmall( 

117 'The cell is not large enough in ' 

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

119 f'Recommended rmax = {recommended_r_max}') 

120 

121 

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

123 recommended_r_max = 5.0 

124 vol = cell.volume 

125 for i in range(3): 

126 if pbc[i]: 

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

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

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

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

131 return recommended_r_max 

132 

133 

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

135 atom2xyz = atoms.get_positions() 

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

137 

138 

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

140 return np.product(get_containing_cell_length(atoms))