Coverage for /builds/ase/ase/ase/geometry/cell.py: 88.17%

93 statements  

« prev     ^ index     » next       coverage.py v7.5.3, created at 2025-08-02 00:12 +0000

1# fmt: off 

2 

3# Copyright (C) 2010, Jesper Friis 

4# (see accompanying license files for details). 

5 

6# XXX bravais objects need to hold tolerance eps, *or* temember variant 

7# from the beginning. 

8# 

9# Should they hold a 'cycle' argument or other data to reconstruct a particular 

10# cell? (E.g. rotation, niggli transform) 

11# 

12# Implement total ordering of Bravais classes 1-14 

13 

14import numpy as np 

15from numpy import arccos, cos, dot, pi, sin, sqrt 

16from numpy.linalg import norm 

17 

18 

19def unit_vector(x): 

20 """Return a unit vector in the same direction as x.""" 

21 y = np.array(x, dtype='float') 

22 return y / norm(y) 

23 

24 

25def angle(x, y): 

26 """Return the angle between vectors a and b in degrees.""" 

27 return arccos(dot(x, y) / (norm(x) * norm(y))) * 180. / pi 

28 

29 

30def cell_to_cellpar(cell, radians=False): 

31 """Returns the cell parameters [a, b, c, alpha, beta, gamma]. 

32 

33 Angles are in degrees unless radian=True is used. 

34 """ 

35 lengths = [np.linalg.norm(v) for v in cell] 

36 angles = [] 

37 for i in range(3): 

38 j = i - 1 

39 k = i - 2 

40 ll = lengths[j] * lengths[k] 

41 if ll > 1e-16: 

42 x = np.dot(cell[j], cell[k]) / ll 

43 angle = 180.0 / pi * arccos(x) 

44 else: 

45 angle = 90.0 

46 angles.append(angle) 

47 if radians: 

48 angles = [angle * pi / 180 for angle in angles] 

49 return np.array(lengths + angles) 

50 

51 

52def cellpar_to_cell(cellpar, ab_normal=(0, 0, 1), a_direction=None): 

53 """Return a 3x3 cell matrix from cellpar=[a,b,c,alpha,beta,gamma]. 

54 

55 Angles must be in degrees. 

56 

57 The returned cell is orientated such that a and b 

58 are normal to `ab_normal` and a is parallel to the projection of 

59 `a_direction` in the a-b plane. 

60 

61 Default `a_direction` is (1,0,0), unless this is parallel to 

62 `ab_normal`, in which case default `a_direction` is (0,0,1). 

63 

64 The returned cell has the vectors va, vb and vc along the rows. The 

65 cell will be oriented such that va and vb are normal to `ab_normal` 

66 and va will be along the projection of `a_direction` onto the a-b 

67 plane. 

68 

69 Example: 

70 

71 >>> from ase.geometry.cell import cellpar_to_cell 

72 

73 >>> cell = cellpar_to_cell([1, 2, 4, 10, 20, 30], (0, 1, 1), (1, 2, 3)) 

74 >>> np.round(cell, 3) 

75 array([[ 0.816, -0.408, 0.408], 

76 [ 1.992, -0.13 , 0.13 ], 

77 [ 3.859, -0.745, 0.745]]) 

78 

79 """ 

80 if a_direction is None: 

81 if np.linalg.norm(np.cross(ab_normal, (1, 0, 0))) < 1e-5: 

82 a_direction = (0, 0, 1) 

83 else: 

84 a_direction = (1, 0, 0) 

85 

86 # Define rotated X,Y,Z-system, with Z along ab_normal and X along 

87 # the projection of a_direction onto the normal plane of Z. 

88 ad = np.array(a_direction) 

89 Z = unit_vector(ab_normal) 

90 X = unit_vector(ad - dot(ad, Z) * Z) 

91 Y = np.cross(Z, X) 

92 

93 # Express va, vb and vc in the X,Y,Z-system 

94 alpha, beta, gamma = 90., 90., 90. 

95 if isinstance(cellpar, (int, float)): 

96 a = b = c = cellpar 

97 elif len(cellpar) == 1: 

98 a = b = c = cellpar[0] 

99 elif len(cellpar) == 3: 

100 a, b, c = cellpar 

101 else: 

102 a, b, c, alpha, beta, gamma = cellpar 

103 

104 # Handle orthorhombic cells separately to avoid rounding errors 

105 eps = 2 * np.spacing(90.0, dtype=np.float64) # around 1.4e-14 

106 # alpha 

107 if abs(abs(alpha) - 90) < eps: 

108 cos_alpha = 0.0 

109 else: 

110 cos_alpha = cos(alpha * pi / 180.0) 

111 # beta 

112 if abs(abs(beta) - 90) < eps: 

113 cos_beta = 0.0 

114 else: 

115 cos_beta = cos(beta * pi / 180.0) 

116 # gamma 

117 if abs(gamma - 90) < eps: 

118 cos_gamma = 0.0 

119 sin_gamma = 1.0 

120 elif abs(gamma + 90) < eps: 

121 cos_gamma = 0.0 

122 sin_gamma = -1.0 

123 else: 

124 cos_gamma = cos(gamma * pi / 180.0) 

125 sin_gamma = sin(gamma * pi / 180.0) 

126 

127 # Build the cell vectors 

128 va = a * np.array([1, 0, 0]) 

129 vb = b * np.array([cos_gamma, sin_gamma, 0]) 

130 cx = cos_beta 

131 cy = (cos_alpha - cos_beta * cos_gamma) / sin_gamma 

132 cz_sqr = 1. - cx * cx - cy * cy 

133 assert cz_sqr >= 0 

134 cz = sqrt(cz_sqr) 

135 vc = c * np.array([cx, cy, cz]) 

136 

137 # Convert to the Cartesian x,y,z-system 

138 abc = np.vstack((va, vb, vc)) 

139 T = np.vstack((X, Y, Z)) 

140 cell = dot(abc, T) 

141 

142 return cell 

143 

144 

145def metric_from_cell(cell): 

146 """Calculates the metric matrix from cell, which is given in the 

147 Cartesian system.""" 

148 cell = np.asarray(cell, dtype=float) 

149 return np.dot(cell, cell.T) 

150 

151 

152def complete_cell(cell): 

153 """Calculate complete cell with missing lattice vectors. 

154 

155 Returns a new 3x3 ndarray. 

156 """ 

157 

158 cell = np.array(cell, dtype=float) 

159 missing = np.nonzero(~cell.any(axis=1))[0] 

160 

161 if len(missing) == 3: 

162 cell.flat[::4] = 1.0 

163 if len(missing) == 2: 

164 # Must decide two vectors: 

165 V, s, WT = np.linalg.svd(cell.T) 

166 sf = [s[0], 1, 1] 

167 cell = (V @ np.diag(sf) @ WT).T 

168 if np.sign(np.linalg.det(cell)) < 0: 

169 cell[missing[0]] = -cell[missing[0]] 

170 elif len(missing) == 1: 

171 i = missing[0] 

172 cell[i] = np.cross(cell[i - 2], cell[i - 1]) 

173 cell[i] /= np.linalg.norm(cell[i]) 

174 

175 return cell 

176 

177 

178def is_orthorhombic(cell): 

179 """Check that cell only has stuff in the diagonal.""" 

180 return not (np.flatnonzero(cell) % 4).any() 

181 

182 

183def orthorhombic(cell): 

184 """Return cell as three box dimensions or raise ValueError.""" 

185 if not is_orthorhombic(cell): 

186 raise ValueError('Not orthorhombic') 

187 return cell.diagonal().copy() 

188 

189 

190# We make the Cell object available for import from here for compatibility 

191from ase.cell import Cell # noqa