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
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3# Copyright (C) 2010, Jesper Friis
4# (see accompanying license files for details).
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
14import numpy as np
15from numpy import arccos, cos, dot, pi, sin, sqrt
16from numpy.linalg import norm
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)
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
30def cell_to_cellpar(cell, radians=False):
31 """Returns the cell parameters [a, b, c, alpha, beta, gamma].
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)
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].
55 Angles must be in degrees.
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.
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).
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.
69 Example:
71 >>> from ase.geometry.cell import cellpar_to_cell
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]])
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)
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)
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
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)
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])
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)
142 return cell
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)
152def complete_cell(cell):
153 """Calculate complete cell with missing lattice vectors.
155 Returns a new 3x3 ndarray.
156 """
158 cell = np.array(cell, dtype=float)
159 missing = np.nonzero(~cell.any(axis=1))[0]
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])
175 return cell
178def is_orthorhombic(cell):
179 """Check that cell only has stuff in the diagonal."""
180 return not (np.flatnonzero(cell) % 4).any()
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()
190# We make the Cell object available for import from here for compatibility
191from ase.cell import Cell # noqa