Coverage for /builds/ase/ase/ase/cell.py: 99.32%
147 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3from typing import Mapping, Sequence, Union
5import numpy as np
7import ase
8from ase.utils import pbc2pbc
9from ase.utils.arraywrapper import arraylike
11__all__ = ['Cell']
14@arraylike
15class Cell:
16 """Parallel epipedal unit cell of up to three dimensions.
18 This object resembles a 3x3 array whose [i, j]-th element is the jth
19 Cartesian coordinate of the ith unit vector.
21 Cells of less than three dimensions are represented by placeholder
22 unit vectors that are zero."""
24 ase_objtype = 'cell' # For JSON'ing
26 def __init__(self, array):
27 """Create cell.
29 Parameters:
31 array: 3x3 arraylike object
32 The three cell vectors: cell[0], cell[1], and cell[2].
33 """
34 array = np.asarray(array, dtype=float)
35 assert array.shape == (3, 3)
36 self.array = array
38 def cellpar(self, radians=False):
39 """Get unit cell parameters. Sequence of 6 numbers.
41 First three are unit cell vector lengths and second three
42 are angles between them::
44 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]
46 in degrees.
48 See also :func:`ase.geometry.cell.cell_to_cellpar`."""
49 from ase.geometry.cell import cell_to_cellpar
50 return cell_to_cellpar(self.array, radians)
52 def todict(self):
53 return dict(array=self.array)
55 @classmethod
56 def ascell(cls, cell):
57 """Return argument as a Cell object. See :meth:`ase.cell.Cell.new`.
59 A new Cell object is created if necessary."""
60 if isinstance(cell, cls):
61 return cell
62 return cls.new(cell)
64 @classmethod
65 def new(cls, cell=None):
66 """Create new cell from any parameters.
68 If cell is three numbers, assume three lengths with right angles.
70 If cell is six numbers, assume three lengths, then three angles.
72 If cell is 3x3, assume three cell vectors."""
74 if cell is None:
75 cell = np.zeros((3, 3))
77 cell = np.array(cell, float)
79 if cell.shape == (3,):
80 cell = np.diag(cell)
81 elif cell.shape == (6,):
82 from ase.geometry.cell import cellpar_to_cell
83 cell = cellpar_to_cell(cell)
84 elif cell.shape != (3, 3):
85 raise ValueError('Cell must be length 3 sequence, length 6 '
86 'sequence or 3x3 matrix!')
88 cellobj = cls(cell)
89 return cellobj
91 @classmethod
92 def fromcellpar(cls, cellpar, ab_normal=(0, 0, 1), a_direction=None):
93 """Return new Cell from cell lengths and angles.
95 See also :func:`~ase.geometry.cell.cellpar_to_cell()`."""
96 from ase.geometry.cell import cellpar_to_cell
97 cell = cellpar_to_cell(cellpar, ab_normal, a_direction)
98 return cls(cell)
100 def get_bravais_lattice(self, eps=2e-4, *, pbc=True):
101 """Return :class:`~ase.lattice.BravaisLattice` for this cell:
103 >>> from ase.cell import Cell
105 >>> cell = Cell.fromcellpar([4, 4, 4, 60, 60, 60])
106 >>> print(cell.get_bravais_lattice())
107 FCC(a=5.65685)
109 .. note:: The Bravais lattice object follows the AFlow
110 conventions. ``cell.get_bravais_lattice().tocell()`` may
111 differ from the original cell by a permutation or other
112 operation which maps it to the AFlow convention. For
113 example, the orthorhombic lattice enforces a < b < c.
115 To build a bandpath for a particular cell, use
116 :meth:`ase.cell.Cell.bandpath` instead of this method.
117 This maps the kpoints back to the original input cell.
119 """
120 from ase.lattice import identify_lattice
121 pbc = self.mask() & pbc2pbc(pbc)
122 lat, _op = identify_lattice(self, eps=eps, pbc=pbc)
123 return lat
125 def bandpath(
126 self,
127 path: str = None,
128 npoints: int = None,
129 *,
130 density: float = None,
131 special_points: Mapping[str, Sequence[float]] = None,
132 eps: float = 2e-4,
133 pbc: Union[bool, Sequence[bool]] = True
134 ) -> "ase.dft.kpoints.BandPath":
135 """Build a :class:`~ase.dft.kpoints.BandPath` for this cell.
137 If special points are None, determine the Bravais lattice of
138 this cell and return a suitable Brillouin zone path with
139 standard special points.
141 If special special points are given, interpolate the path
142 directly from the available data.
144 Parameters:
146 path: string
147 String of special point names defining the path, e.g. 'GXL'.
148 npoints: int
149 Number of points in total. Note that at least one point
150 is added for each special point in the path.
151 density: float
152 density of kpoints along the path in Å⁻¹.
153 special_points: dict
154 Dictionary mapping special points to scaled kpoint coordinates.
155 For example ``{'G': [0, 0, 0], 'X': [1, 0, 0]}``.
156 eps: float
157 Tolerance for determining Bravais lattice.
158 pbc: three bools
159 Whether cell is periodic in each direction. Normally not
160 necessary. If cell has three nonzero cell vectors, use
161 e.g. pbc=[1, 1, 0] to request a 2D bandpath nevertheless.
163 Example
164 -------
166 >>> from ase.cell import Cell
168 >>> cell = Cell.fromcellpar([4, 4, 4, 60, 60, 60])
169 >>> cell.bandpath('GXW', npoints=20)
170 BandPath(path='GXW', cell=[3x3], special_points={GKLUWX}, kpts=[20x3])
172 """
173 # TODO: Combine with the rotation transformation from bandpath()
175 cell = self.uncomplete(pbc)
177 if special_points is None:
178 from ase.lattice import identify_lattice
179 lat, op = identify_lattice(cell, eps=eps)
180 bandpath = lat.bandpath(path, npoints=npoints, density=density)
181 return bandpath.transform(op)
182 else:
183 from ase.dft.kpoints import BandPath, resolve_custom_points
184 path, special_points = resolve_custom_points(
185 path, special_points, eps=eps)
186 bandpath = BandPath(cell, path=path, special_points=special_points)
187 return bandpath.interpolate(npoints=npoints, density=density)
189 def uncomplete(self, pbc):
190 """Return new cell, zeroing cell vectors where not periodic."""
191 _pbc = np.empty(3, bool)
192 _pbc[:] = pbc
193 cell = self.copy()
194 cell[~_pbc] = 0
195 return cell
197 def complete(self):
198 """Convert missing cell vectors into orthogonal unit vectors."""
199 from ase.geometry.cell import complete_cell
200 return Cell(complete_cell(self.array))
202 def copy(self):
203 """Return a copy of this cell."""
204 return Cell(self.array.copy())
206 def mask(self):
207 """Boolean mask of which cell vectors are nonzero."""
208 return self.any(1)
210 @property
211 def rank(self) -> int:
212 """"Return the dimension of the cell.
214 Equal to the number of nonzero lattice vectors."""
215 # The name ndim clashes with ndarray.ndim
216 return sum(self.mask())
218 @property
219 def orthorhombic(self) -> bool:
220 """Return whether this cell is represented by a diagonal matrix."""
221 from ase.geometry.cell import is_orthorhombic
222 return is_orthorhombic(self)
224 def lengths(self):
225 """Return the length of each lattice vector as an array."""
226 return np.linalg.norm(self, axis=1)
228 def angles(self):
229 """Return an array with the three angles alpha, beta, and gamma."""
230 return self.cellpar()[3:].copy()
232 def __array__(self, dtype=None, copy=False):
233 return np.array(self.array, dtype=dtype, copy=copy)
235 def __bool__(self):
236 return bool(self.any()) # need to convert from np.bool_
238 @property
239 def volume(self) -> float:
240 """Get the volume of this cell.
242 If there are less than 3 lattice vectors, return 0."""
243 # Fail or 0 for <3D cells?
244 # Definitely 0 since this is currently a property.
245 # I think normally it is more convenient just to get zero
246 return np.abs(np.linalg.det(self))
248 @property
249 def handedness(self) -> int:
250 """Sign of the determinant of the matrix of cell vectors.
252 1 for right-handed cells, -1 for left, and 0 for cells that
253 do not span three dimensions."""
254 return int(np.sign(np.linalg.det(self)))
256 def scaled_positions(self, positions) -> np.ndarray:
257 """Calculate scaled positions from Cartesian positions.
259 The scaled positions are the positions given in the basis
260 of the cell vectors. For the purpose of defining the basis, cell
261 vectors that are zero will be replaced by unit vectors as per
262 :meth:`~ase.cell.Cell.complete`."""
263 return np.linalg.solve(self.complete().T, np.transpose(positions)).T
265 def cartesian_positions(self, scaled_positions) -> np.ndarray:
266 """Calculate Cartesian positions from scaled positions."""
267 return scaled_positions @ self.complete()
269 def reciprocal(self) -> 'Cell':
270 """Get reciprocal lattice as a Cell object.
272 The reciprocal cell is defined such that
274 cell.reciprocal() @ cell.T == np.diag(cell.mask())
276 within machine precision.
278 Does not include factor of 2 pi."""
279 icell = Cell(np.linalg.pinv(self).transpose())
280 icell[~self.mask()] = 0.0 # type: ignore[index]
281 return icell
283 def normal(self, i):
284 """Normal vector of the two vectors with index different from i.
286 This is the cross product of those vectors in cyclic order from i."""
287 return np.cross(self[i - 2], self[i - 1])
289 def normals(self):
290 """Normal vectors of each axis as a 3x3 matrix."""
291 return np.array([self.normal(i) for i in range(3)])
293 def area(self, i):
294 """Area spanned by the two vectors with index different from i."""
295 return np.linalg.norm(self.normal(i))
297 def areas(self):
298 """Areas spanned by cell vector pairs (1, 2), (2, 0), and (0, 2)."""
299 return np.linalg.norm(self.normals(), axis=1)
301 def __repr__(self):
302 if self.orthorhombic:
303 numbers = self.lengths().tolist()
304 else:
305 numbers = self.tolist()
307 return f'Cell({numbers})'
309 def niggli_reduce(self, eps=1e-5):
310 """Niggli reduce this cell, returning a new cell and mapping.
312 See also :func:`ase.build.tools.niggli_reduce_cell`."""
313 from ase.build.tools import niggli_reduce_cell
314 cell, op = niggli_reduce_cell(self, epsfactor=eps)
315 result = Cell(cell)
316 return result, op
318 def minkowski_reduce(self):
319 """Minkowski-reduce this cell, returning new cell and mapping.
321 See also :func:`ase.geometry.minkowski_reduction.minkowski_reduce`."""
322 from ase.geometry.minkowski_reduction import minkowski_reduce
323 cell, op = minkowski_reduce(self, self.mask())
324 result = Cell(cell)
325 return result, op
327 def permute_axes(self, permutation):
328 """Permute axes of cell."""
329 assert (np.sort(permutation) == np.arange(3)).all()
330 permuted = Cell(self[permutation][:, permutation])
331 return permuted
333 def standard_form(self, form='lower'):
334 """Rotate axes such that unit cell is lower/upper triangular. The cell
335 handedness is preserved.
337 A lower-triangular cell with positive diagonal entries is a canonical
338 (i.e. unique) description. For a left-handed cell the diagonal entries
339 are negative.
341 Parameters:
343 form: str
344 'lower' or 'upper' triangular form. The default is 'lower'.
346 Returns:
348 rcell: the standardized cell object
350 Q: ndarray
351 The orthogonal transformation. Here, rcell @ Q = cell, where cell
352 is the input cell and rcell is the lower triangular (output) cell.
353 """
355 sign = self.handedness
356 if sign == 0:
357 sign = 1
359 # LQ decomposition provides an axis-aligned description of the cell.
360 # Q is an orthogonal matrix and L is a lower triangular matrix. The
361 # decomposition is a unique description if the diagonal elements are
362 # all positive (negative for a left-handed cell).
363 if form == 'lower':
364 Q, L = np.linalg.qr(self.T)
365 Q = Q.T
366 L = L.T
367 elif form == 'upper':
368 Q, L = np.linalg.qr(self.T[::-1, ::-1])
369 Q = Q.T[::-1, ::-1]
370 L = L.T[::-1, ::-1]
371 else:
372 raise ValueError('form must be either "lower" or "upper"')
374 # correct the signs of the diagonal elements
375 signs = np.sign(np.diag(L))
376 indices = np.where(signs == 0)[0]
377 signs[indices] = 1
378 indices = np.where(signs != sign)[0]
379 L[:, indices] *= -1
380 Q[indices] *= -1
381 return Cell(L), Q
383 # XXX We want a reduction function that brings the cell into
384 # standard form as defined by Setyawan and Curtarolo.