Coverage for dynasor / qpoints / tools.py: 98%

1import itertools 

2from fractions import Fraction 

3 

4import numpy as np 

5from numpy.typing import NDArray 

6 

7from dynasor.modes.tools import inv 

8 

9 

10def get_supercell_qpoints_along_path( 

11 path: list[tuple[str, str]], 

12 coordinates: dict[str, NDArray[float]], 

13 primitive_cell: NDArray[float], 

14 super_cell: NDArray[float]) -> list[NDArray[float]]: 

15 r""" 

16 Returns the q-points commensurate with the given supercell along the specific path. 

17 

18 Parameters 

19 ---------- 

20 path 

21 List of pairs of q-point labels. 

22 coordinates 

23 Dict with q-point labels and coordinates as keys and values, respectively; 

24 there must be one entry for each q-point label used in :attr:`path`. 

25 primitive_cell 

26 Cell metric of the primitive cell with lattice vectors as rows. 

27 super_cell 

28 Cell metric of the supercell with lattice vectors as rows. 

29 

30 Returns 

31 ------- 

32 A list of the accessible q-point coordinates along the specified segment. 

33 

34 Example 

35 -------- 

36 The following example illustrates how to retrieve the q-points that 

37 can be sampled using a supercell comprising :math:`6 \times 6 \times 6` 

38 conventional (4-atom) unit cells of FCC Al along the path X-:math:`\Gamma`-L. 

39 

40 >>> import numpy as np 

41 >>> from ase.build import bulk 

42 >>> from dynasor.qpoints import get_supercell_qpoints_along_path 

43 >>> prim = bulk('Al', 'fcc', a=4.0) 

44 >>> supercell = bulk('Al', 'fcc', a=4.0, cubic=True).repeat(6) 

45 >>> path = [('X', 'G'), ('G', 'L'), ('L', 'W')] 

46 >>> coordinates = dict(X=[0.5, 0.5, 0], G=[0, 0, 0], 

47 ... L=[0.5, 0.5, 0.5], W=[0.5, 0.25, 0.75]) 

48 >>> qpoints = get_supercell_qpoints_along_path( 

49 ... path, coordinates, prim.cell, supercell.cell) 

50 

51 """ 

52 from .lattice import Lattice 

53 lat = Lattice(primitive_cell, super_cell) 

54 

55 for lbl in np.array(path).flatten(): 

56 if lbl not in coordinates: 56 ↛ 57line 56 didn't jump to line 57 because the condition on line 56 was never true

57 raise ValueError(f'Unknown point in path: {lbl}') 

58 

59 # build the segments 

60 supercell_paths = [] 

61 for k, (l1, l2) in enumerate(path): 

62 q1 = np.array(coordinates[l1], dtype=float) 

63 q2 = np.array(coordinates[l2], dtype=float) 

64 dynasor_path, _ = lat.make_path(q1, q2) 

65 supercell_paths.append(dynasor_path) 

66 return supercell_paths 

67 

68 

69def find_on_line( 

70 start: NDArray[float], 

71 stop: NDArray[float], 

72 P: NDArray[int], 

73) -> set[NDArray[float]]: 

74 """Find fractional distances between start and stop combatible with P. 

75 

76 A supercell is defined by `P @ c = S` for some repetition matrix `P` and we 

77 want to find fractions so that:: 

78 

79 [start + f * (stop - start)] @ P = n 

80 

81 Parameters 

82 ---------- 

83 start 

84 Start of line in reduced supercell coordinates. 

85 stop 

86 End of line in reduced supercell coordinates. 

87 P 

88 Repetion matrix defining the supercell. 

89 """ 

90 

91 if np.allclose(start, stop): 

92 return [Fraction(0, 1)] 

93 

94 start = np.array([Fraction(s).limit_denominator() for s in start]) 

95 stop = np.array([Fraction(s).limit_denominator() for s in stop]) 

96 

97 A = start @ P 

98 B = (stop - start) @ P 

99 

100 fracs = None 

101 for a, b in zip(A, B): 

102 fs = solve_Diophantine(a, b) 

103 if fs is None: # "inf" solutions 

104 continue 

105 elif fs == []: # No solutions 

106 return [] 

107 fracs = set(fs) if fracs is None else fracs.intersection(fs) 

108 return sorted(fracs) 

109 

110 

111def solve_Diophantine(a: Fraction, b: Fraction) -> list[Fraction]: 

112 """Solve n = a + xb for all n in Z and a,b in Q such that 0 <= x <= 1.""" 

113 

114 if b == 0: 

115 if a.denominator == 1: 

116 return None 

117 else: 

118 return [] 

119 

120 if b < 0: 

121 right = np.ceil(a) 

122 left = np.floor(a + b) 

123 else: 

124 left = np.floor(a) 

125 right = np.ceil(a + b) 

126 

127 ns = np.arange(left, right + 1) 

128 fracs = [Fraction(n - a, b) for n in ns] 

129 fracs = [f for f in fracs if 0 <= f <= 1] 

130 

131 return fracs 

132 

133 

134def get_commensurate_lattice_points(P: NDArray[int]) -> NDArray[float]: 

135 """Return commensurate points for a supercell defined by repetition matrix `P`. 

136 

137 Finds all `n` such that `n = f P` where `f` is between 0 and 1 

138 

139 Parameters 

140 ---------- 

141 P 

142 The repetion matrix relating the primitive and supercell. 

143 

144 Returns 

145 ------- 

146 The commensurate lattice points. 

147 """ 

148 inv_P_matrix = inv(P, as_fraction=True) 

149 

150 assert np.all(P @ inv_P_matrix == np.eye(3)) 

151 

152 n_max = np.where(P > 0, P, 0).sum(axis=0) + 1 

153 n_min = np.where(P < 0, P, 0).sum(axis=0) 

154 

155 ranges = [np.arange(*n) for n in zip(n_min, n_max)] 

156 

157 lattice_points = [] 

158 for lp in itertools.product(*ranges): 

159 tmp = lp @ inv_P_matrix 

160 if np.all(tmp >= 0) and np.all(tmp < 1): 

161 lattice_points.append(lp) 

162 

163 assert len(lattice_points) == len(set(lattice_points)) 

164 lattice_points = np.array(lattice_points) 

165 return lattice_points