Coverage for local_installation/dynasor/tools/structures.py: 96%
47 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-08-05 09:53 +0000
« prev ^ index » next coverage.py v7.3.2, created at 2024-08-05 09:53 +0000
1import numpy as np
2from ase import Atoms
5def align_structure(atoms: Atoms, atol: float = 1e-5):
6 """
8 Rotate and realign atoms object such that
9 * the first cell vector points along the x-directon
10 * the second cell vector lies in the xy-plane
12 Modifies input ``atoms`` object in place.
14 Parameters
15 ----------
16 atoms
17 input structure to be rotated aligned wtih the x,y,z coordinte system
18 atol
19 absolute tolerance used for sanity checking the cell
20 """
21 _align_a_onto_xy(atoms, atol)
22 _align_a_onto_x(atoms, atol)
23 _align_b_onto_xy(atoms, atol)
26def _align_a_onto_xy(atoms, atol):
27 """Rotate cell so that a is in the xy-plane"""
29 # get angle towards xy
30 # will break if a is along z
31 assert np.any(atoms.cell[0, :2])
33 cell = atoms.cell.array.copy()
35 a = cell[0]
36 a_xy = a.copy()
37 a_xy[2] = 0 # projection of a onto xy-plane
39 # angle between a and xy-plane
40 cosa = np.dot(a, a_xy) / np.linalg.norm(a) / np.linalg.norm(a_xy)
42 # cosa should be in the interval (0, 1]
43 assert not np.isclose(cosa, 0)
44 if cosa > 1: 44 ↛ 45line 44 didn't jump to line 45, because the condition on line 44 was never true
45 assert np.isclose(cosa, 1)
46 cosa = min(cosa, 1)
47 cosa = max(cosa, 0)
49 # angle between a and xy-plane in degs
50 angle_xy_deg = np.rad2deg(np.arccos(cosa))
52 # get unit vector to rotate around
53 vec = np.cross(a_xy, [0, 0, 1])
54 vec = vec / np.linalg.norm(vec)
55 assert vec[2] == 0
57 # Determine if the rotation should be positive or negative depending on
58 # whether a is pointing in the +z or -z direction
59 sign = -1 if a[2] > 0 else +1
61 # rotate
62 atoms.rotate(sign * angle_xy_deg, vec, rotate_cell=True)
64 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0), atoms.cell
67def _align_a_onto_x(atoms, atol):
68 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0) # make sure a is in xy-plane
70 a = atoms.cell[0]
71 a_x = a[0]
72 a_y = a[1]
74 # angle between a and x-axis (a is already in xy-plane)
76 # tan = y / x -> angle = arctan y / x "=" atan2(y, x)
77 angle_rad = np.arctan2(a_y, a_x)
78 angle_deg = np.rad2deg(angle_rad)
80 atoms.rotate(-angle_deg, [0, 0, 1], rotate_cell=True)
82 assert np.isclose(atoms.cell[0, 1], 0, atol=atol, rtol=0), atoms.cell
83 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0), atoms.cell
86def _align_b_onto_xy(atoms, atol):
87 assert np.isclose(atoms.cell[0, 1], 0, atol=atol, rtol=0) # make sure a is along x
88 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0) # make sure a is along x
90 # rotate so that b is in xy plane
91 # project b onto the yz-plane
92 b = atoms.cell[1]
93 b_y = b[1]
94 b_z = b[2]
95 angle_rad = np.arctan2(b_z, b_y)
96 angle_deg = np.rad2deg(angle_rad)
98 atoms.rotate(-angle_deg, [1, 0, 0], rotate_cell=True)
100 assert np.isclose(atoms.cell[0, 1], 0, atol=atol, rtol=0) # make sure a is in xy-plane
101 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0) # make sure a is in xy-plane
102 assert np.isclose(atoms.cell[1, 2], 0, atol=atol, rtol=0), atoms.cell