Coverage for ase / cell.py: 99.32%
147 statements
« prev ^ index » next coverage.py v7.13.3, created at 2026-02-04 10:20 +0000
1"""Cell."""
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
30 ----------
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
51 return cell_to_cellpar(self.array, radians)
53 def todict(self):
54 return dict(array=self.array)
56 @classmethod
57 def ascell(cls, cell):
58 """Return argument as a Cell object. See :meth:`ase.cell.Cell.new`.
60 A new Cell object is created if necessary."""
61 if isinstance(cell, cls):
62 return cell
63 return cls.new(cell)
65 @classmethod
66 def new(cls, cell=None):
67 """Create new cell from any parameters.
69 If cell is three numbers, assume three lengths with right angles.
71 If cell is six numbers, assume three lengths, then three angles.
73 If cell is 3x3, assume three cell vectors."""
75 if cell is None:
76 cell = np.zeros((3, 3))
78 cell = np.array(cell, float)
80 if cell.shape == (3,):
81 cell = np.diag(cell)
82 elif cell.shape == (6,):
83 from ase.geometry.cell import cellpar_to_cell
85 cell = cellpar_to_cell(cell)
86 elif cell.shape != (3, 3):
87 raise ValueError(
88 'Cell must be length 3 sequence, length 6 '
89 'sequence or 3x3 matrix!'
90 )
92 cellobj = cls(cell)
93 return cellobj
95 @classmethod
96 def fromcellpar(cls, cellpar, ab_normal=(0, 0, 1), a_direction=None):
97 """Return new Cell from cell lengths and angles.
99 See also :func:`~ase.geometry.cell.cellpar_to_cell()`."""
100 from ase.geometry.cell import cellpar_to_cell
102 cell = cellpar_to_cell(cellpar, ab_normal, a_direction)
103 return cls(cell)
105 def get_bravais_lattice(self, eps=2e-4, *, pbc=True):
106 """Return :class:`~ase.lattice.BravaisLattice` for this cell:
108 >>> from ase.cell import Cell
110 >>> cell = Cell.fromcellpar([4, 4, 4, 60, 60, 60])
111 >>> print(cell.get_bravais_lattice())
112 FCC(a=5.65685)
114 .. note:: The Bravais lattice object follows the AFlow
115 conventions. ``cell.get_bravais_lattice().tocell()`` may
116 differ from the original cell by a permutation or other
117 operation which maps it to the AFlow convention. For
118 example, the orthorhombic lattice enforces a < b < c.
120 To build a bandpath for a particular cell, use
121 :meth:`ase.cell.Cell.bandpath` instead of this method.
122 This maps the kpoints back to the original input cell.
124 """
125 from ase.lattice import identify_lattice
127 pbc = self.mask() & pbc2pbc(pbc)
128 lat, _op = identify_lattice(self, eps=eps, pbc=pbc)
129 return lat
131 def bandpath(
132 self,
133 path: str = None,
134 npoints: int = None,
135 *,
136 density: float = None,
137 special_points: Mapping[str, Sequence[float]] = None,
138 eps: float = 2e-4,
139 pbc: Union[bool, Sequence[bool]] = True,
140 ) -> 'ase.dft.kpoints.BandPath':
141 """Build a :class:`~ase.dft.kpoints.BandPath` for this cell.
143 If special points are None, determine the Bravais lattice of
144 this cell and return a suitable Brillouin zone path with
145 standard special points.
147 If special special points are given, interpolate the path
148 directly from the available data.
150 Parameters
151 ----------
152 path: string
153 String of special point names defining the path, e.g. 'GXL'.
154 npoints: int
155 Number of points in total. Note that at least one point
156 is added for each special point in the path.
157 density: float
158 density of kpoints along the path in Å⁻¹.
159 special_points: dict
160 Dictionary mapping special points to scaled kpoint coordinates.
161 For example ``{'G': [0, 0, 0], 'X': [1, 0, 0]}``.
162 eps: float
163 Tolerance for determining Bravais lattice.
164 pbc: three bools
165 Whether cell is periodic in each direction. Normally not
166 necessary. If cell has three nonzero cell vectors, use
167 e.g. pbc=[1, 1, 0] to request a 2D bandpath nevertheless.
169 Example
170 -------
172 >>> from ase.cell import Cell
174 >>> cell = Cell.fromcellpar([4, 4, 4, 60, 60, 60])
175 >>> cell.bandpath('GXW', npoints=20)
176 BandPath(path='GXW', cell=[3x3], special_points={GKLUWX}, kpts=[20x3])
178 """
179 # TODO: Combine with the rotation transformation from bandpath()
181 cell = self.uncomplete(pbc)
183 if special_points is None:
184 from ase.lattice import identify_lattice
186 lat, op = identify_lattice(cell, eps=eps)
187 bandpath = lat.bandpath(path, npoints=npoints, density=density)
188 return bandpath.transform(op)
189 else:
190 from ase.dft.kpoints import BandPath, resolve_custom_points
192 path, special_points = resolve_custom_points(
193 path, special_points, eps=eps
194 )
195 bandpath = BandPath(cell, path=path, special_points=special_points)
196 return bandpath.interpolate(npoints=npoints, density=density)
198 def uncomplete(self, pbc):
199 """Return new cell, zeroing cell vectors where not periodic."""
200 _pbc = np.empty(3, bool)
201 _pbc[:] = pbc
202 cell = self.copy()
203 cell[~_pbc] = 0
204 return cell
206 def complete(self):
207 """Convert missing cell vectors into orthogonal unit vectors."""
208 from ase.geometry.cell import complete_cell
210 return Cell(complete_cell(self.array))
212 def copy(self):
213 """Return a copy of this cell."""
214 return Cell(self.array.copy())
216 def mask(self):
217 """Boolean mask of which cell vectors are nonzero."""
218 return self.any(1)
220 @property
221 def rank(self) -> int:
222 """Return the dimension of the cell.
224 Equal to the number of nonzero lattice vectors."""
225 # The name ndim clashes with ndarray.ndim
226 return sum(self.mask())
228 @property
229 def orthorhombic(self) -> bool:
230 """Return whether this cell is represented by a diagonal matrix."""
231 from ase.geometry.cell import is_orthorhombic
233 return is_orthorhombic(self)
235 def lengths(self):
236 """Return the length of each lattice vector as an array."""
237 return np.linalg.norm(self, axis=1)
239 def angles(self):
240 """Return an array with the three angles alpha, beta, and gamma."""
241 return self.cellpar()[3:].copy()
243 def __array__(self, dtype=None, copy=False):
244 return np.array(self.array, dtype=dtype, copy=copy)
246 def __bool__(self):
247 return bool(self.any()) # need to convert from np.bool_
249 @property
250 def volume(self) -> float:
251 """Get the volume of this cell.
253 If there are less than 3 lattice vectors, return 0."""
254 # Fail or 0 for <3D cells?
255 # Definitely 0 since this is currently a property.
256 # I think normally it is more convenient just to get zero
257 return np.abs(np.linalg.det(self))
259 @property
260 def handedness(self) -> int:
261 """Sign of the determinant of the matrix of cell vectors.
263 1 for right-handed cells, -1 for left, and 0 for cells that
264 do not span three dimensions."""
265 return int(np.sign(np.linalg.det(self)))
267 def scaled_positions(self, positions) -> np.ndarray:
268 """Calculate scaled positions from Cartesian positions.
270 The scaled positions are the positions given in the basis
271 of the cell vectors. For the purpose of defining the basis, cell
272 vectors that are zero will be replaced by unit vectors as per
273 :meth:`~ase.cell.Cell.complete`."""
274 return np.linalg.solve(self.complete().T, np.transpose(positions)).T
276 def cartesian_positions(self, scaled_positions) -> np.ndarray:
277 """Calculate Cartesian positions from scaled positions."""
278 return scaled_positions @ self.complete()
280 def reciprocal(self) -> 'Cell':
281 """Get reciprocal lattice as a Cell object.
283 The reciprocal cell is defined such that
285 cell.reciprocal() @ cell.T == np.diag(cell.mask())
287 within machine precision.
289 Does not include factor of 2 pi."""
290 icell = Cell(np.linalg.pinv(self).transpose())
291 icell[~self.mask()] = 0.0 # type: ignore[index]
292 return icell
294 def normal(self, i):
295 """Normal vector of the two vectors with index different from i.
297 This is the cross product of those vectors in cyclic order from i."""
298 return np.cross(self[i - 2], self[i - 1])
300 def normals(self):
301 """Normal vectors of each axis as a 3x3 matrix."""
302 return np.array([self.normal(i) for i in range(3)])
304 def area(self, i):
305 """Area spanned by the two vectors with index different from i."""
306 return np.linalg.norm(self.normal(i))
308 def areas(self):
309 """Areas spanned by cell vector pairs (1, 2), (2, 0), and (0, 1)."""
310 return np.linalg.norm(self.normals(), axis=1)
312 def __repr__(self):
313 if self.orthorhombic:
314 numbers = self.lengths().tolist()
315 else:
316 numbers = self.tolist()
318 return f'Cell({numbers})'
320 def niggli_reduce(self, eps=1e-5):
321 """Niggli reduce this cell, returning a new cell and mapping.
323 See also :func:`ase.build.tools.niggli_reduce_cell`."""
324 from ase.build.tools import niggli_reduce_cell
326 cell, op = niggli_reduce_cell(self, epsfactor=eps)
327 result = Cell(cell)
328 return result, op
330 def minkowski_reduce(self):
331 """Minkowski-reduce this cell, returning new cell and mapping.
333 See also :func:`ase.geometry.minkowski_reduction.minkowski_reduce`."""
334 from ase.geometry.minkowski_reduction import minkowski_reduce
336 cell, op = minkowski_reduce(self, self.mask())
337 result = Cell(cell)
338 return result, op
340 def permute_axes(self, permutation):
341 """Permute axes of cell."""
342 assert (np.sort(permutation) == np.arange(3)).all()
343 permuted = Cell(self[permutation][:, permutation])
344 return permuted
346 def standard_form(self, form='lower'):
347 """Rotate axes such that unit cell is lower/upper triangular. The cell
348 handedness is preserved.
350 A lower-triangular cell with positive diagonal entries is a canonical
351 (i.e. unique) description. For a left-handed cell the diagonal entries
352 are negative.
354 Parameters
355 ----------
356 form: str
357 'lower' or 'upper' triangular form. The default is 'lower'.
359 Returns
360 -------
361 rcell: the standardized cell object
363 Q: ndarray
364 The orthogonal transformation. Here, rcell @ Q = cell, where cell
365 is the input cell and rcell is the lower triangular (output) cell.
366 """
368 sign = self.handedness
369 if sign == 0:
370 sign = 1
372 # LQ decomposition provides an axis-aligned description of the cell.
373 # Q is an orthogonal matrix and L is a lower triangular matrix. The
374 # decomposition is a unique description if the diagonal elements are
375 # all positive (negative for a left-handed cell).
376 if form == 'lower':
377 Q, L = np.linalg.qr(self.T)
378 Q = Q.T
379 L = L.T
380 elif form == 'upper':
381 Q, L = np.linalg.qr(self.T[::-1, ::-1])
382 Q = Q.T[::-1, ::-1]
383 L = L.T[::-1, ::-1]
384 else:
385 raise ValueError('form must be either "lower" or "upper"')
387 # correct the signs of the diagonal elements
388 signs = np.sign(np.diag(L))
389 indices = np.where(signs == 0)[0]
390 signs[indices] = 1
391 indices = np.where(signs != sign)[0]
392 L[:, indices] *= -1
393 Q[indices] *= -1
394 return Cell(L), Q
396 # XXX We want a reduction function that brings the cell into
397 # standard form as defined by Setyawan and Curtarolo.