Coverage for local_installation/dynasor/qpoints/tools.py: 99%

1import itertools 

2from fractions import Fraction 

3from typing import Dict, List, Tuple 

4 

5import numpy as np 

6from numpy.typing import NDArray 

7from ase import Atoms 

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 supercell_paths 

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

34 

35 Example 

36 -------- 

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

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

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

40 

41 >>> import numpy as np 

42 >>> from ase.build import bulk 

43 >>> from dynasor.qpoints import get_supercell_qpoints_along_path 

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

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

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

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

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

49 >>> qpoints = get_supercell_qpoints_along_path(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(start: NDArray, stop: NDArray, P: NDArray): 

70 """Find fractional distances between start and stop combatible with P 

71 

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

73 want to find fractions so that 

74 

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

76 

77 Parameters 

78 ---------- 

79 start 

80 start of line in reduced supercell coordinates 

81 stop 

82 end of line in reduced supercell coordinates 

83 P 

84 repetion matrix defining the supercell 

85 """ 

86 

87 if np.allclose(start, stop): 

88 return [Fraction(0, 1)] 

89 

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

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

92 

93 A = start @ P 

94 B = (stop - start) @ P 

95 

96 fracs = None 

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

98 fs = solve_Diophantine(a, b) 

99 if fs is None: # "inf" solutions 

100 continue 

101 elif fs == []: # No solutions 

102 return [] 

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

104 return sorted(fracs) 

105 

106 

107def solve_Diophantine(a: Fraction, b: Fraction) -> List[Fraction]: 

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

109 

110 if b == 0: 

111 if a.denominator == 1: 

112 return None 

113 else: 

114 return [] 

115 

116 if b < 0: 

117 right = np.ceil(a) 

118 left = np.floor(a + b) 

119 else: 

120 left = np.floor(a) 

121 right = np.ceil(a + b) 

122 

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

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

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

126 

127 return fracs 

128 

129 

130def det(A): 

131 """Determinant of an integer matrix using Laplace cofactor expansion""" 

132 if len(A) == 2: 

133 return A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0] 

134 d = 0 

135 for i, B in enumerate(A[0]): # along first row 

136 minor = np.hstack([A[1:, :i], A[1:, i+1:]]) 

137 d += (-1)**i * B * det(minor) 

138 assert np.isclose(d, np.linalg.det(A)) 

139 return d 

140 

141 

142def inv(A): 

143 """Takes the inverse of an integer 3x3 matrix based on Cayley-Hamilton""" 

144 

145 detx2 = det(A) * 2 # Denominator "determinant times two" 

146 

147 # Numerator 

148 numerator = ((np.trace(A)**2 - np.trace(A @ A)) * np.diag([1, 1, 1]) 

149 - 2 * A * np.trace(A) 

150 + 2 * A @ A) 

151 

152 # We want the sign to be in the Numerator 

153 if detx2 < 0: 

154 detx2 = -detx2 

155 numerator = -numerator 

156 

157 inverse = numerator / detx2 

158 assert np.allclose(inverse, np.linalg.inv(A)) 

159 

160 # Return inverse, numerator (int matrix) and denominator (int) 

161 return inverse, numerator, detx2 

162 

163 

164def get_P_matrix(c, S): 

165 """ P c = S -> c.T P.T = S.T 

166 

167 The P matrix must be an integer matrix 

168 """ 

169 PT = np.linalg.solve(c.T, S.T) 

170 P_float = PT.T 

171 P = np.round(P_float).astype(int) 

172 if not np.allclose(P_float, P) or not np.allclose(P @ c, S): 

173 raise ValueError( 

174 f'Please check that the supercell metric ({S}) is related to the' 

175 f' the primitive cell {c} by an integer transformation matrix.') 

176 return P 

177 

178 

179def get_commensurate_lattice_points(P: NDArray) -> NDArray: 

180 """Return commensurate points for a supercell defined by repetition matrix P 

181 

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

183 

184 Parameters 

185 ---------- 

186 P 

187 the repetion matrix relating the primitive and supercell 

188 

189 Returns 

190 ------- 

191 lattice_points 

192 the commensurate lattice points 

193 """ 

194 

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

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

197 

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

199 

200 inv_P_matrix, num, den = inv(P) 

201 

202 lattice_points = [] 

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

204 s = lp @ num # here we skip the denominator to keep everything integer 

205 # the denominator is also integer so no numerics here 

206 if np.all(s >= 0) and np.all(s < den): 

207 lattice_points.append(lp) 

208 

209 lattice_points = np.array(lattice_points) 

210 

211 # Begin sane checks... 

212 

213 # No duplicates 

214 assert len(lattice_points) == len(np.unique(lattice_points, axis=0)) 

215 

216 # Did we get everyone? 

217 assert len(lattice_points) == abs(det(P)) 

218 

219 return lattice_points 

220 

221 

222def get_index_offset(supercell: Atoms, prim: Atoms, atol=1e-3, rtol=0.0): 

223 """ 

224 Get the basis index and primitive cell offsets for a supercell 

225 """ 

226 index, offset = [], [] 

227 for pos in supercell.positions: 

228 spos = np.linalg.solve(prim.cell.T, pos) 

229 for i, spos2 in enumerate(prim.get_scaled_positions()): 

230 off = spos - spos2 

231 off_round = np.round(off) 

232 if not np.allclose(off, off_round, atol=atol, rtol=rtol): 

233 continue 

234 index.append(i) 

235 off = off_round.astype(int) 

236 assert np.allclose(off, off_round) 

237 offset.append(off) 

238 break 

239 else: 

240 raise ValueError('prim not compatible with atoms') 

241 

242 index, offset = np.array(index), np.array(offset) 

243 return index, offset