r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import math

2from typing import List, Optional, Tuple, Union

4import numpy as np

5from ase import Atoms

6from ase.cell import Cell

9class CellTooSmall(Exception):

10 pass

13class VolumeNotDefined(Exception):

14 pass

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.

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.

31 Parameters:

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.

38 nbins : int

39 Number of bins to divide the rdf into.

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.

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.

50 no_dists : bool

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

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

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

63 check_cell_and_r_max(atoms, rmax)

65 dm = distance_matrix

66 if dm is None:

67 dm = atoms.get_all_distances(mic=True)

69 rdf = np.zeros(nbins + 1)

70 dr = float(rmax / nbins)

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

73 natoms = len(atoms)

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)

80 indices_triu = np.triu(indices)

81 for index in range(nbins + 1):

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

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

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

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

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

98 if no_dists:

99 return rdf[1:]

101 return rdf[1:], rr

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

105 cell = atoms.get_cell()

106 pbc = atoms.get_pbc()

108 vol = atoms.cell.volume

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

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

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

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

140 return np.product(get_containing_cell_length(atoms))