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