Coverage for /builds/kinetik161/ase/ase/cell.py: 99.30%
143 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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=float):
231 if dtype != float:
232 raise ValueError('Cannot convert cell to array of type {}'
233 .format(dtype))
234 return self.array
236 def __bool__(self):
237 return bool(self.any()) # need to convert from np.bool_
239 @property
240 def volume(self) -> float:
241 """Get the volume of this cell.
243 If there are less than 3 lattice vectors, return 0."""
244 # Fail or 0 for <3D cells?
245 # Definitely 0 since this is currently a property.
246 # I think normally it is more convenient just to get zero
247 return np.abs(np.linalg.det(self))
249 @property
250 def handedness(self) -> int:
251 """Sign of the determinant of the matrix of cell vectors.
253 1 for right-handed cells, -1 for left, and 0 for cells that
254 do not span three dimensions."""
255 return int(np.sign(np.linalg.det(self)))
257 def scaled_positions(self, positions) -> np.ndarray:
258 """Calculate scaled positions from Cartesian positions.
260 The scaled positions are the positions given in the basis
261 of the cell vectors. For the purpose of defining the basis, cell
262 vectors that are zero will be replaced by unit vectors as per
263 :meth:`~ase.cell.Cell.complete`."""
264 return np.linalg.solve(self.complete().T, np.transpose(positions)).T
266 def cartesian_positions(self, scaled_positions) -> np.ndarray:
267 """Calculate Cartesian positions from scaled positions."""
268 return scaled_positions @ self.complete()
270 def reciprocal(self) -> 'Cell':
271 """Get reciprocal lattice as a Cell object.
273 The reciprocal cell is defined such that
275 cell.reciprocal() @ cell.T == np.diag(cell.mask())
277 within machine precision.
279 Does not include factor of 2 pi."""
280 icell = Cell(np.linalg.pinv(self).transpose())
281 icell[~self.mask()] = 0.0 # type: ignore[index]
282 return icell
284 def normal(self, i):
285 """Normal vector of the two vectors with index different from i.
287 This is the cross product of those vectors in cyclic order from i."""
288 return np.cross(self[i - 2], self[i - 1])
290 def normals(self):
291 """Normal vectors of each axis as a 3x3 matrix."""
292 return np.array([self.normal(i) for i in range(3)])
294 def area(self, i):
295 """Area spanned by the two vectors with index different from i."""
296 return np.linalg.norm(self.normal(i))
298 def areas(self):
299 """Areas spanned by cell vector pairs (1, 2), (2, 0), and (0, 2)."""
300 return np.linalg.norm(self.normals(), axis=1)
302 def __repr__(self):
303 if self.orthorhombic:
304 numbers = self.lengths().tolist()
305 else:
306 numbers = self.tolist()
308 return f'Cell({numbers})'
310 def niggli_reduce(self, eps=1e-5):
311 """Niggli reduce this cell, returning a new cell and mapping.
313 See also :func:`ase.build.tools.niggli_reduce_cell`."""
314 from ase.build.tools import niggli_reduce_cell
315 cell, op = niggli_reduce_cell(self, epsfactor=eps)
316 result = Cell(cell)
317 return result, op
319 def minkowski_reduce(self):
320 """Minkowski-reduce this cell, returning new cell and mapping.
322 See also :func:`ase.geometry.minkowski_reduction.minkowski_reduce`."""
323 from ase.geometry.minkowski_reduction import minkowski_reduce
324 cell, op = minkowski_reduce(self, self.mask())
325 result = Cell(cell)
326 return result, op
328 def permute_axes(self, permutation):
329 """Permute axes of cell."""
330 assert (np.sort(permutation) == np.arange(3)).all()
331 permuted = Cell(self[permutation][:, permutation])
332 return permuted
334 def standard_form(self):
335 """Rotate axes such that unit cell is lower triangular. The cell
336 handedness is preserved.
338 A lower-triangular cell with positive diagonal entries is a canonical
339 (i.e. unique) description. For a left-handed cell the diagonal entries
340 are negative.
342 Returns:
344 rcell: the standardized cell object
346 Q: ndarray
347 The orthogonal transformation. Here, rcell @ Q = cell, where cell
348 is the input cell and rcell is the lower triangular (output) cell.
349 """
351 sign = self.handedness
352 if sign == 0:
353 sign = 1
355 # LQ decomposition provides an axis-aligned description of the cell.
356 # Q is an orthogonal matrix and L is a lower triangular matrix. The
357 # decomposition is a unique description if the diagonal elements are
358 # all positive (negative for a left-handed cell).
359 Q, L = np.linalg.qr(self.T)
360 Q = Q.T
361 L = L.T
363 # correct the signs of the diagonal elements
364 signs = np.sign(np.diag(L))
365 indices = np.where(signs == 0)[0]
366 signs[indices] = 1
367 indices = np.where(signs != sign)[0]
368 L[:, indices] *= -1
369 Q[indices] *= -1
370 return Cell(L), Q
372 # XXX We want a reduction function that brings the cell into
373 # standard form as defined by Setyawan and Curtarolo.