Coverage for dynasor/qpoints/tools.py: 98%
67 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-03 19:46 +0000
1import itertools
2from fractions import Fraction
4import numpy as np
5from numpy.typing import NDArray
7from dynasor.modes.tools import inv
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 A list of the accessible q-point coordinates along the specified segment.
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.
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)
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(
70 start: NDArray[float],
71 stop: NDArray[float],
72 P: NDArray[int],
73) -> list[Fraction]:
74 """Find fractional distances between start and stop compatible with P.
76 A supercell is defined by `P @ c = S` for some repetition matrix `P` and we
77 want to find fractions so that::
79 [start + f * (stop - start)] @ P = n
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 Repetition matrix defining the supercell.
89 """
91 if np.allclose(start, stop):
92 return [Fraction(0, 1)]
94 start = np.array([Fraction(s).limit_denominator() for s in start])
95 stop = np.array([Fraction(s).limit_denominator() for s in stop])
97 A = start @ P
98 B = (stop - start) @ P
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)
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."""
114 if b == 0:
115 if a.denominator == 1:
116 return None
117 else:
118 return []
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)
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]
131 return fracs
134def is_qpoint_commensurate(
135 q_point: NDArray[float],
136 cell: NDArray[float],
137 atol: float = 1e-5,
138) -> bool:
139 r"""Return whether a single q-point is commensurate with a periodic cell.
141 A q-point is commensurate with the cell when it is a reciprocal-lattice
142 vector of that cell, i.e. when its reduced coordinates
143 :math:`m_i = \mathbf{q}\cdot\mathbf{a}_i / 2\pi` are all integers, where the
144 :math:`\mathbf{a}_i` are the cell vectors (rows of :attr:`cell`). Only then
145 is :math:`\rho(\mathbf{q}) = \sum_j e^{i\mathbf{q}\cdot\mathbf{r}_j}`
146 independent of which periodic image is chosen for each atom; at other
147 q-points the correlation functions can contain artifacts.
149 The reduced coordinates are formed with the full cell matrix, so the test is
150 correct for non-orthogonal cells (a per-axis check would not be).
152 Parameters
153 ----------
154 q_point
155 A single q-point ``(3,)`` in Cartesian coordinates (rad/Å).
156 cell
157 Cell metric with lattice vectors as rows (Å).
158 atol
159 Absolute tolerance on the reduced coordinates when testing for integers.
160 """
161 q_point = np.asarray(q_point, dtype=float)
162 cell = np.asarray(cell, dtype=float)
163 reduced = q_point @ cell.T / (2 * np.pi) # m_i = q . a_i / 2pi
164 return bool(np.all(np.abs(reduced - np.rint(reduced)) < atol))
167def get_commensurate_lattice_points(P: NDArray[int]) -> NDArray[float]:
168 """Return commensurate points for a supercell defined by repetition matrix `P`.
170 Finds all `n` such that `n = f P` where `f` is between 0 and 1
172 Parameters
173 ----------
174 P
175 The repetition matrix relating the primitive and supercell.
177 Returns
178 -------
179 The commensurate lattice points.
180 """
181 inv_P_matrix = inv(P, as_fraction=True)
183 assert np.all(P @ inv_P_matrix == np.eye(3))
185 n_max = np.where(P > 0, P, 0).sum(axis=0) + 1
186 n_min = np.where(P < 0, P, 0).sum(axis=0)
188 ranges = [np.arange(*n) for n in zip(n_min, n_max)]
190 lattice_points = []
191 for lp in itertools.product(*ranges):
192 tmp = lp @ inv_P_matrix
193 if np.all(tmp >= 0) and np.all(tmp < 1):
194 lattice_points.append(lp)
196 assert len(lattice_points) == len(set(lattice_points))
197 lattice_points = np.array(lattice_points)
198 return lattice_points