# Coverage for local_installation/dynasor/tools/structures.py: 96%

## 47 statements

, 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

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)

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]