 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# Copyright (C) 2010, Jesper Friis

2# (see accompanying license files for details).

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

5# from the beginning.

6#

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

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

9#

10# Implement total ordering of Bravais classes 1-14

12import numpy as np

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

14from numpy.linalg import norm

17def unit_vector(x):

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

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

20 return y / norm(y)

23def angle(x, y):

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

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

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

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

32 """

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

34 angles = []

35 for i in range(3):

36 j = i - 1

37 k = i - 2

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

39 if ll > 1e-16:

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

41 angle = 180.0 / pi * arccos(x)

42 else:

43 angle = 90.0

44 angles.append(angle)

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

47 return np.array(lengths + angles)

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

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

53 Angles must be in degrees.

55 The returned cell is orientated such that a and b

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

57 `a_direction` in the a-b plane.

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

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

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

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

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

65 plane.

67 Example:

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

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

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

72 [ 1.992, -0.13 , 0.13 ],

73 [ 3.859, -0.745, 0.745]])

75 """

76 if a_direction is None:

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

78 a_direction = (0, 0, 1)

79 else:

80 a_direction = (1, 0, 0)

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

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

85 Z = unit_vector(ab_normal)

87 Y = np.cross(Z, X)

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

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

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

92 a = b = c = cellpar

93 elif len(cellpar) == 1:

94 a = b = c = cellpar

95 elif len(cellpar) == 3:

96 a, b, c = cellpar

97 else:

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

100 # Handle orthorhombic cells separately to avoid rounding errors

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

102 # alpha

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

104 cos_alpha = 0.0

105 else:

106 cos_alpha = cos(alpha * pi / 180.0)

107 # beta

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

109 cos_beta = 0.0

110 else:

111 cos_beta = cos(beta * pi / 180.0)

112 # gamma

113 if abs(gamma - 90) < eps:

114 cos_gamma = 0.0

115 sin_gamma = 1.0

116 elif abs(gamma + 90) < eps:

117 cos_gamma = 0.0

118 sin_gamma = -1.0

119 else:

120 cos_gamma = cos(gamma * pi / 180.0)

121 sin_gamma = sin(gamma * pi / 180.0)

123 # Build the cell vectors

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

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

126 cx = cos_beta

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

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

129 assert cz_sqr >= 0

130 cz = sqrt(cz_sqr)

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

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

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

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

136 cell = dot(abc, T)

138 return cell

141def metric_from_cell(cell):

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

143 Cartesian system."""

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

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

148def complete_cell(cell):

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

151 Returns a new 3x3 ndarray.

152 """

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

155 missing = np.nonzero(~cell.any(axis=1))

157 if len(missing) == 3:

158 cell.flat[::4] = 1.0

159 if len(missing) == 2:

160 # Must decide two vectors:

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

162 sf = [s, 1, 1]

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

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

165 cell[missing] = -cell[missing]

166 elif len(missing) == 1:

167 i = missing

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

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

171 return cell

174def is_orthorhombic(cell):

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

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

179def orthorhombic(cell):

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

181 if not is_orthorhombic(cell):

182 raise ValueError('Not orthorhombic')

183 return cell.diagonal().copy()

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

187from ase.cell import Cell # noqa