Coverage for local_installation/dynasor/qpoints/tools.py: 97%
107 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-12-21 12:02 +0000
1import itertools
2from fractions import Fraction
3from typing import Dict, List, Tuple
5import numpy as np
6from numpy.typing import NDArray
7from ase import Atoms
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.
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
30 Returns
31 -------
32 supercell_paths
33 A list of the accessible q-point coordinates along the specified segment
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.
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)
51 """
52 from .lattice import Lattice
53 lat = Lattice(primitive_cell, super_cell)
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}')
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
69def find_on_line(start: NDArray, stop: NDArray, P: NDArray):
70 """Find fractional distances between start and stop combatible with P
72 A supercell is defined by P @ c = S for some repetition matrix P and we
73 want to find fractions so that
75 [start + f * (stop - start)] @ P = n
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 """
87 if np.allclose(start, stop):
88 return [Fraction(0, 1)]
90 start = np.array([Fraction(s).limit_denominator() for s in start])
91 stop = np.array([Fraction(s).limit_denominator() for s in stop])
93 A = start @ P
94 B = (stop - start) @ P
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)
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"""
110 if b == 0:
111 if a.denominator == 1:
112 return None
113 else:
114 return []
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)
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]
127 return fracs
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
142def inv(A):
143 """Takes the inverse of an integer 3x3 matrix based on Cayley-Hamilton"""
145 detx2 = det(A) * 2 # Denominator "determinant times two"
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)
152 # We want the sign to be in the Numerator
153 if detx2 < 0:
154 detx2 = -detx2
155 numerator = -numerator
157 inverse = numerator / detx2
158 assert np.allclose(inverse, np.linalg.inv(A))
160 # Return inverse, numerator (int matrix) and denominator (int)
161 return inverse, numerator, detx2
164def get_P_matrix(c, S):
165 """ P c = S -> c.T P.T = S.T
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
179def get_commensurate_lattice_points(P: NDArray) -> NDArray:
180 """Return commensurate points for a supercell defined by repetition matrix P
182 Finds all n such that n = f P where f is between 0 and 1
184 Parameters
185 ----------
186 P
187 the repetion matrix relating the primitive and supercell
189 Returns
190 -------
191 lattice_points
192 the commensurate lattice points
193 """
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)
198 ranges = [np.arange(*n) for n in zip(n_min, n_max)]
200 inv_P_matrix, num, den = inv(P)
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)
209 lattice_points = np.array(lattice_points)
211 # Begin sane checks...
213 # No duplicates
214 assert len(lattice_points) == len(np.unique(lattice_points, axis=0))
216 # Did we get everyone?
217 assert len(lattice_points) == abs(det(P))
219 return lattice_points
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 """
227 if len(prim) > len(supercell): 227 ↛ 228line 227 didn't jump to line 228, because the condition on line 227 was never true
228 raise ValueError('prim contains more atoms than supercell')
230 index, offset = [], []
231 for pos in supercell.positions:
232 spos = np.linalg.solve(prim.cell.T, pos)
233 for i, spos2 in enumerate(prim.get_scaled_positions()):
234 off = spos - spos2
235 off_round = np.round(off)
236 if not np.allclose(off, off_round, atol=atol, rtol=rtol):
237 continue
238 index.append(i)
239 off = off_round.astype(int)
240 assert np.allclose(off, off_round)
241 offset.append(off)
242 break
243 else:
244 raise ValueError('prim not compatible with supercell')
246 index, offset = np.array(index), np.array(offset)
247 return index, offset