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

1import numpy as np 

2from ase import Atoms 

3 

4 

5def align_structure(atoms: Atoms, atol: float = 1e-5): 

6 """ 

7 

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 

11 

12 Modifies input ``atoms`` object in place. 

13 

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) 

24 

25 

26def _align_a_onto_xy(atoms, atol): 

27 """Rotate cell so that a is in the xy-plane""" 

28 

29 # get angle towards xy 

30 # will break if a is along z 

31 assert np.any(atoms.cell[0, :2]) 

32 

33 cell = atoms.cell.array.copy() 

34 

35 a = cell[0] 

36 a_xy = a.copy() 

37 a_xy[2] = 0 # projection of a onto xy-plane 

38 

39 # angle between a and xy-plane 

40 cosa = np.dot(a, a_xy) / np.linalg.norm(a) / np.linalg.norm(a_xy) 

41 

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) 

48 

49 # angle between a and xy-plane in degs 

50 angle_xy_deg = np.rad2deg(np.arccos(cosa)) 

51 

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 

56 

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 

60 

61 # rotate 

62 atoms.rotate(sign * angle_xy_deg, vec, rotate_cell=True) 

63 

64 assert np.isclose(atoms.cell[0, 2], 0, atol=atol, rtol=0), atoms.cell 

65 

66 

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 

69 

70 a = atoms.cell[0] 

71 a_x = a[0] 

72 a_y = a[1] 

73 

74 # angle between a and x-axis (a is already in xy-plane) 

75 

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) 

79 

80 atoms.rotate(-angle_deg, [0, 0, 1], rotate_cell=True) 

81 

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 

84 

85 

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 

89 

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) 

97 

98 atoms.rotate(-angle_deg, [1, 0, 0], rotate_cell=True) 

99 

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