Coverage for dynasor / tools / structures.py: 88%
111 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-16 12:31 +0000
1from typing import Optional
2import numpy as np
3from ase import Atoms
4from ase.geometry import get_distances
5from ase.geometry import find_mic
6from dynasor.modes.atoms import Prim
7from numpy.typing import NDArray
10def get_displacements(atoms: Atoms,
11 atoms_ideal: Atoms,
12 check_mic: Optional[bool] = True,
13 cell_tol: Optional[float] = 1e-4) -> NDArray[float]:
14 """Returns the the smallest possible displacements between a
15 displaced configuration relative to an ideal (reference)
16 configuration.
18 Parameters
19 ----------
20 atoms
21 Structure with displaced atoms.
22 ideal
23 Ideal configuration relative to which displacements are computed.
24 check_mic
25 Whether to check minimum image convention.
26 cell_tol
27 Cell tolerance; if cell missmatch more than tol value error is raised.
28 """
30 if not np.array_equal(atoms.numbers, atoms_ideal.numbers): 30 ↛ 31line 30 didn't jump to line 31 because the condition on line 30 was never true
31 raise ValueError('Atomic numbers do not match.')
32 if np.linalg.norm(atoms.cell - atoms_ideal.cell) > cell_tol: 32 ↛ 33line 32 didn't jump to line 33 because the condition on line 32 was never true
33 raise ValueError('Cells do not match.')
35 u = atoms.positions - atoms_ideal.positions
36 return get_displacements_from_u(u, atoms_ideal.cell, check_mic=True)
39def get_displacements_from_u(
40 u: NDArray[float],
41 cell: NDArray[float],
42 check_mic: Optional[bool] = True,
43) -> NDArray[float]:
44 """wraps displacements using mic"""
45 if check_mic: 45 ↛ 47line 45 didn't jump to line 47 because the condition on line 45 was always true
46 u, _ = find_mic(u, cell)
47 return u
50def find_permutation(atoms: Atoms, atoms_ref: Atoms) -> list[int]:
51 """ Returns the best permutation of atoms for mapping one
52 configuration onto another.
54 Parameters
55 ----------
56 atoms
57 Configuration to be permuted.
58 atoms_ref
59 Configuration onto which to map.
61 Example
62 -------
63 After obtaining the permutation via ``p = find_permutation(atoms1, atoms2)``
64 the reordered structure ``atoms1[p]`` will give the closest match
65 to ``atoms2``.
66 """
67 if np.linalg.norm(atoms.cell - atoms_ref.cell) > 1e-4: 67 ↛ 68line 67 didn't jump to line 68 because the condition on line 67 was never true
68 raise ValueError('Cells do not match.')
70 permutation = []
71 for i in range(len(atoms_ref)):
72 dist_row = get_distances(
73 atoms.positions, atoms_ref.positions[i], cell=atoms_ref.cell, pbc=True)[1][:, 0]
74 permutation.append(np.argmin(dist_row))
76 if len(set(permutation)) != len(permutation): 76 ↛ 77line 76 didn't jump to line 77 because the condition on line 76 was never true
77 raise Exception('Duplicates in permutation.')
78 for i, p in enumerate(permutation):
79 if atoms[p].symbol != atoms_ref[i].symbol: 79 ↛ 80line 79 didn't jump to line 80 because the condition on line 79 was never true
80 raise Exception('Matching lattice sites have different occupation.')
81 return permutation
84def align_structure(atoms: Atoms, atol: Optional[float] = 1e-5) -> None:
85 """
86 Rotates and realigns a structure such that
87 * the first cell vector points along the x-directon
88 * the second cell vector lies in the xy-plane
90 Note that this function modifies the :attr:`atoms` object in place.
92 Parameters
93 ----------
94 atoms
95 Input structure to be rotated aligned with the x,y,z coordinte system.
96 atol
97 Absolute tolerance used for sanity checking the cell.
98 """
99 _align_a_onto_xy(atoms, atol)
100 _align_a_onto_x(atoms, atol)
101 _align_b_onto_xy(atoms, atol)
104def get_offset_index(
105 primitive: Atoms,
106 supercell: Atoms,
107 tol: Optional[float] = 0.01,
108 wrap: Optional[bool] = True,
109) -> tuple[NDArray[float], NDArray[float]]:
110 """ Returns the basis index and primitive cell offsets for a supercell.
112 This implementation uses a simple iteration procedure that should be fairly quick.
113 If more stability is needed consider the following approach:
115 * find the P-matrix: `P = ideal.cell @ prim.cell_inv.T`
116 * compensate for strain: `P *= len(ideal)/len(prim)/det(P)`
117 * generate the reference structure: `ref_atoms = make_supercell(round(P), prim)`
118 * find the assignment using `ref_atoms` via the Hungarian algorithm using the mic distances
120 Parameters
121 ----------
122 primitive
123 Primitive cell.
124 supercell
125 Some ideal repetition of the primitive cell.
126 tol
127 Tolerance length parameter. Increase to allow for slgihtly rattled or strained cells.
128 wrap
129 It might happen that the ideal cell boundary cuts through a unit cell
130 whose lattice points lie inside the ideal cell. If there is a basis, an
131 atom belonging to this unit cell might get wrapped while another is
132 not. Then the wrapped atom now belongs to a lattice point outside the P
133 matrix so to say. This would result in more lattice points than
134 expected from `N_unit = len(ideal)/len(prim)`.
136 Returns
137 -------
138 offsets
139 The lattice points as integers in `(N, 3)`-array.
140 index
141 The basis indices as integers in `(N,)`-array.
142 """
144 if not isinstance(primitive, Atoms): 144 ↛ 145line 144 didn't jump to line 145 because the condition on line 144 was never true
145 raise ValueError
146 if not isinstance(supercell, Atoms): 146 ↛ 147line 146 didn't jump to line 147 because the condition on line 146 was never true
147 raise ValueError
149 prim = Prim(primitive)
151 from dynasor.modes.tools import inv
153 P = get_P_matrix(primitive.cell, supercell.cell) # P C = S
154 P_inv = inv(P)
156 lattice, basis = [], []
157 # Pick an atom in the supercell
158 for pos_ideal in supercell.positions:
159 # Does this atom perhaps belong to site "index"?
160 for index, pos_prim in enumerate(primitive.positions):
161 # if so we can remove the basis position vector and should end up on a lattice site
162 diff_pos = pos_ideal - pos_prim
163 # The lattice site has integer coordinates in reduced coordinates
164 prim_spos = diff_pos @ prim.inv_cell.T
165 # Rounding should not affect the coordinate much if it is integer
166 prim_spos_round = np.round(prim_spos).astype(int)
167 # If the rounded spos and unrounded spos are the same
168 if np.allclose(prim_spos, prim_spos_round, rtol=0, atol=tol):
169 # Since P_inv is represented using fractions we can neatly
170 # write the supercell spos of the lattice point using fractions
171 # and easily determine if it needs wrapping or not without
172 # worry about numerics
173 ideal_spos = prim_spos_round @ P_inv
174 # wrap if needed
175 ideal_spos_wrap = ideal_spos % 1 if wrap else ideal_spos
176 # This should be integer again
177 prim_spos_wrap = (ideal_spos_wrap @ P).astype(int)
178 # add the results and break out from the basis site loop
179 lattice.append(prim_spos_wrap)
180 basis.append(index)
181 break
182 else: # we get here by not breaking out from the basis site loop.
183 # This means that the candidate lattice site where not close to integers
185 raise ValueError(f' {prim_spos} {prim_spos_round} Supercell not compatible '
186 'with primitive cell.')
188 lattice = np.array(lattice)
189 basis = np.array(basis)
191 # We should have found len(ideal) unique positions
192 lattice_basis = [tuple((*lp, i)) for lp, i in zip(lattice, basis)]
193 assert len(set(lattice_basis)) == len(supercell)
195 return lattice, basis
198def get_P_matrix(
199 c: NDArray[float],
200 S: NDArray[float],
201) -> NDArray[float]:
202 """Returns the P matrix, i.e., the `3x3` integer matrix :math:`P` that satisfies
204 .. math::
206 P c = S
208 Here, :math:`c` is the primitive cell metric and :math:`S` is the
209 supercell metric as row vectors. Note that the above condition is
210 equivalent to:
212 .. math::
214 c^T P^T = S^T
216 Parameters
217 ----------
218 c
219 Cell metric of the primitive structure.
220 S
221 Cell metric of the supercell.
222 """
223 PT = np.linalg.solve(c.T, S.T)
224 P_float = PT.T
225 P = np.round(P_float).astype(int)
226 if not np.allclose(P_float, P) or not np.allclose(P @ c, S):
227 raise ValueError(
228 f'Please check that the supercell metric ({S}) is related to the'
229 f' the primitive cell {c} by an integer transformation matrix.')
230 return P
233def _align_a_onto_xy(atoms: Atoms, atol: float) -> None:
234 """ Rotate cell so that a is in the xy-plane. """
236 # get angle towards xy
237 # will break if a is along z
238 assert np.any(atoms.cell[0, :2])
240 cell = atoms.cell.array.copy()
242 a = cell[0]
243 a_xy = a.copy()
244 a_xy[2] = 0 # projection of a onto xy-plane
246 # angle between a and xy-plane
247 cosa = np.dot(a, a_xy) / np.linalg.norm(a) / np.linalg.norm(a_xy)
249 # cosa should be in the interval (0, 1]
250 assert not np.isclose(cosa, 0)
251 if cosa > 1: 251 ↛ 252line 251 didn't jump to line 252 because the condition on line 251 was never true
252 assert np.isclose(cosa, 1)
253 cosa = min(cosa, 1)
254 cosa = max(cosa, 0)
256 # angle between a and xy-plane in degs
257 angle_xy_deg = np.rad2deg(np.arccos(cosa))
259 # get unit vector to rotate around
260 vec = np.cross(a_xy, [0, 0, 1])
261 vec = vec / np.linalg.norm(vec)
262 assert vec[2] == 0
264 # Determine if the rotation should be positive or negative depending on
265 # whether a is pointing in the +z or -z direction
266 sign = -1 if a[2] > 0 else +1
268 # rotate
269 atoms.rotate(sign * angle_xy_deg, vec, rotate_cell=True)
271 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0), atoms.cell
274def _align_a_onto_x(atoms: Atoms, atol: float) -> None:
275 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0) # make sure a is in xy-plane
277 a = atoms.cell[0]
278 a_x = a[0]
279 a_y = a[1]
281 # angle between a and x-axis (a is already in xy-plane)
283 # tan = y / x -> angle = arctan y / x "=" atan2(y, x)
284 angle_rad = np.arctan2(a_y, a_x)
285 angle_deg = np.rad2deg(angle_rad)
287 atoms.rotate(-angle_deg, [0, 0, 1], rotate_cell=True)
289 assert np.isclose(atoms.cell[0, 1], 0, atol=atol, rtol=0), atoms.cell
290 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0), atoms.cell
293def _align_b_onto_xy(atoms: Atoms, atol: float) -> None:
294 assert np.isclose(atoms.cell[0, 1], 0, atol=atol, rtol=0) # make sure a is along x
295 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0) # make sure a is along x
297 # rotate so that b is in xy plane
298 # project b onto the yz-plane
299 b = atoms.cell[1]
300 b_y = b[1]
301 b_z = b[2]
302 angle_rad = np.arctan2(b_z, b_y)
303 angle_deg = np.rad2deg(angle_rad)
305 atoms.rotate(-angle_deg, [1, 0, 0], rotate_cell=True)
307 assert np.isclose(atoms.cell[0, 1], 0, atol=atol, rtol=0) # make sure a is in xy-plane
308 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0) # make sure a is in xy-plane
309 assert np.isclose(atoms.cell[1, 2], 0, atol=atol, rtol=0), atoms.cell