Coverage for /builds/ase/ase/ase/cell.py : 99.30%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import ase
2from typing import Mapping, Sequence, Union
3import numpy as np
4from ase.utils.arraywrapper import arraylike
5from ase.utils import pbc2pbc
8__all__ = ['Cell']
11@arraylike
12class Cell:
13 """Parallel epipedal unit cell of up to three dimensions.
15 This object resembles a 3x3 array whose [i, j]-th element is the jth
16 Cartesian coordinate of the ith unit vector.
18 Cells of less than three dimensions are represented by placeholder
19 unit vectors that are zero."""
21 ase_objtype = 'cell' # For JSON'ing
23 def __init__(self, array):
24 """Create cell.
26 Parameters:
28 array: 3x3 arraylike object
29 The three cell vectors: cell[0], cell[1], and cell[2].
30 """
31 array = np.asarray(array, dtype=float)
32 assert array.shape == (3, 3)
33 self.array = array
35 def cellpar(self, radians=False):
36 """Get unit cell parameters. Sequence of 6 numbers.
38 First three are unit cell vector lengths and second three
39 are angles between them::
41 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]
43 in degrees.
45 See also :func:`ase.geometry.cell.cell_to_cellpar`."""
46 from ase.geometry.cell import cell_to_cellpar
47 return cell_to_cellpar(self.array, radians)
49 def todict(self):
50 return dict(array=self.array)
52 @classmethod
53 def ascell(cls, cell):
54 """Return argument as a Cell object. See :meth:`ase.cell.Cell.new`.
56 A new Cell object is created if necessary."""
57 if isinstance(cell, cls):
58 return cell
59 return cls.new(cell)
61 @classmethod
62 def new(cls, cell=None):
63 """Create new cell from any parameters.
65 If cell is three numbers, assume three lengths with right angles.
67 If cell is six numbers, assume three lengths, then three angles.
69 If cell is 3x3, assume three cell vectors."""
71 if cell is None:
72 cell = np.zeros((3, 3))
74 cell = np.array(cell, float)
76 if cell.shape == (3,):
77 cell = np.diag(cell)
78 elif cell.shape == (6,):
79 from ase.geometry.cell import cellpar_to_cell
80 cell = cellpar_to_cell(cell)
81 elif cell.shape != (3, 3):
82 raise ValueError('Cell must be length 3 sequence, length 6 '
83 'sequence or 3x3 matrix!')
85 cellobj = cls(cell)
86 return cellobj
88 @classmethod
89 def fromcellpar(cls, cellpar, ab_normal=(0, 0, 1), a_direction=None):
90 """Return new Cell from cell lengths and angles.
92 See also :func:`~ase.geometry.cell.cellpar_to_cell()`."""
93 from ase.geometry.cell import cellpar_to_cell
94 cell = cellpar_to_cell(cellpar, ab_normal, a_direction)
95 return cls(cell)
97 def get_bravais_lattice(self, eps=2e-4, *, pbc=True):
98 """Return :class:`~ase.lattice.BravaisLattice` for this cell:
100 >>> cell = Cell.fromcellpar([4, 4, 4, 60, 60, 60])
101 >>> print(cell.get_bravais_lattice())
102 FCC(a=5.65685)
104 .. note:: The Bravais lattice object follows the AFlow
105 conventions. ``cell.get_bravais_lattice().tocell()`` may
106 differ from the original cell by a permutation or other
107 operation which maps it to the AFlow convention. For
108 example, the orthorhombic lattice enforces a < b < c.
110 To build a bandpath for a particular cell, use
111 :meth:`ase.cell.Cell.bandpath` instead of this method.
112 This maps the kpoints back to the original input cell.
114 """
115 from ase.lattice import identify_lattice
116 pbc = self.mask() & pbc2pbc(pbc)
117 lat, op = identify_lattice(self, eps=eps, pbc=pbc)
118 return lat
120 def bandpath(
121 self,
122 path: str = None,
123 npoints: int = None,
124 *,
125 density: float = None,
126 special_points: Mapping[str, Sequence[float]] = None,
127 eps: float = 2e-4,
128 pbc: Union[bool, Sequence[bool]] = True
129 ) -> "ase.dft.kpoints.BandPath":
130 """Build a :class:`~ase.dft.kpoints.BandPath` for this cell.
132 If special points are None, determine the Bravais lattice of
133 this cell and return a suitable Brillouin zone path with
134 standard special points.
136 If special special points are given, interpolate the path
137 directly from the available data.
139 Parameters:
141 path: string
142 String of special point names defining the path, e.g. 'GXL'.
143 npoints: int
144 Number of points in total. Note that at least one point
145 is added for each special point in the path.
146 density: float
147 density of kpoints along the path in Å⁻¹.
148 special_points: dict
149 Dictionary mapping special points to scaled kpoint coordinates.
150 For example ``{'G': [0, 0, 0], 'X': [1, 0, 0]}``.
151 eps: float
152 Tolerance for determining Bravais lattice.
153 pbc: three bools
154 Whether cell is periodic in each direction. Normally not
155 necessary. If cell has three nonzero cell vectors, use
156 e.g. pbc=[1, 1, 0] to request a 2D bandpath nevertheless.
158 Example
159 -------
160 >>> cell = Cell.fromcellpar([4, 4, 4, 60, 60, 60])
161 >>> cell.bandpath('GXW', npoints=20)
162 BandPath(path='GXW', cell=[3x3], special_points={GKLUWX}, kpts=[20x3])
164 """
165 # TODO: Combine with the rotation transformation from bandpath()
167 cell = self.uncomplete(pbc)
169 if special_points is None:
170 from ase.lattice import identify_lattice
171 lat, op = identify_lattice(cell, eps=eps)
172 bandpath = lat.bandpath(path, npoints=npoints, density=density)
173 return bandpath.transform(op)
174 else:
175 from ase.dft.kpoints import BandPath, resolve_custom_points
176 path, special_points = resolve_custom_points(
177 path, special_points, eps=eps)
178 bandpath = BandPath(cell, path=path, special_points=special_points)
179 return bandpath.interpolate(npoints=npoints, density=density)
181 def uncomplete(self, pbc):
182 """Return new cell, zeroing cell vectors where not periodic."""
183 _pbc = np.empty(3, bool)
184 _pbc[:] = pbc
185 cell = self.copy()
186 cell[~_pbc] = 0
187 return cell
189 def complete(self):
190 """Convert missing cell vectors into orthogonal unit vectors."""
191 from ase.geometry.cell import complete_cell
192 return Cell(complete_cell(self.array))
194 def copy(self):
195 """Return a copy of this cell."""
196 return Cell(self.array.copy())
198 def mask(self):
199 """Boolean mask of which cell vectors are nonzero."""
200 return self.any(1)
202 @property
203 def rank(self) -> int:
204 """"Return the dimension of the cell.
206 Equal to the number of nonzero lattice vectors."""
207 # The name ndim clashes with ndarray.ndim
208 return sum(self.mask()) # type: ignore
210 @property
211 def orthorhombic(self) -> bool:
212 """Return whether this cell is represented by a diagonal matrix."""
213 from ase.geometry.cell import is_orthorhombic
214 return is_orthorhombic(self)
216 def lengths(self):
217 """Return the length of each lattice vector as an array."""
218 return np.linalg.norm(self, axis=1)
220 def angles(self):
221 """Return an array with the three angles alpha, beta, and gamma."""
222 return self.cellpar()[3:].copy()
224 def __array__(self, dtype=float):
225 if dtype != float:
226 raise ValueError('Cannot convert cell to array of type {}'
227 .format(dtype))
228 return self.array
230 def __bool__(self):
231 return bool(self.any()) # need to convert from np.bool_
233 @property
234 def volume(self) -> float:
235 """Get the volume of this cell.
237 If there are less than 3 lattice vectors, return 0."""
238 # Fail or 0 for <3D cells?
239 # Definitely 0 since this is currently a property.
240 # I think normally it is more convenient just to get zero
241 return np.abs(np.linalg.det(self))
243 @property
244 def handedness(self) -> int:
245 """Sign of the determinant of the matrix of cell vectors.
247 1 for right-handed cells, -1 for left, and 0 for cells that
248 do not span three dimensions."""
249 return int(np.sign(np.linalg.det(self)))
251 def scaled_positions(self, positions) -> np.ndarray:
252 """Calculate scaled positions from Cartesian positions.
254 The scaled positions are the positions given in the basis
255 of the cell vectors. For the purpose of defining the basis, cell
256 vectors that are zero will be replaced by unit vectors as per
257 :meth:`~ase.cell.Cell.complete`."""
258 return np.linalg.solve(self.complete().T, np.transpose(positions)).T
260 def cartesian_positions(self, scaled_positions) -> np.ndarray:
261 """Calculate Cartesian positions from scaled positions."""
262 return scaled_positions @ self.complete()
264 def reciprocal(self) -> 'Cell':
265 """Get reciprocal lattice as a Cell object.
267 The reciprocal cell is defined such that
269 cell.reciprocal() @ cell.T == np.diag(cell.mask())
271 within machine precision.
273 Does not include factor of 2 pi."""
274 icell = Cell(np.linalg.pinv(self).transpose())
275 icell[~self.mask()] = 0.0 # type: ignore
276 return icell
278 def normal(self, i):
279 """Normal vector of the two vectors with index different from i.
281 This is the cross product of those vectors in cyclic order from i."""
282 return np.cross(self[i - 2], self[i - 1])
284 def normals(self):
285 """Normal vectors of each axis as a 3x3 matrix."""
286 return np.array([self.normal(i) for i in range(3)])
288 def area(self, i):
289 """Area spanned by the two vectors with index different from i."""
290 return np.linalg.norm(self.normal(i))
292 def areas(self):
293 """Areas spanned by cell vector pairs (1, 2), (2, 0), and (0, 2)."""
294 return np.linalg.norm(self.normals(), axis=1)
296 def __repr__(self):
297 if self.orthorhombic:
298 numbers = self.lengths().tolist()
299 else:
300 numbers = self.tolist()
302 return 'Cell({})'.format(numbers)
304 def niggli_reduce(self, eps=1e-5):
305 """Niggli reduce this cell, returning a new cell and mapping.
307 See also :func:`ase.build.tools.niggli_reduce_cell`."""
308 from ase.build.tools import niggli_reduce_cell
309 cell, op = niggli_reduce_cell(self, epsfactor=eps)
310 result = Cell(cell)
311 return result, op
313 def minkowski_reduce(self):
314 """Minkowski-reduce this cell, returning new cell and mapping.
316 See also :func:`ase.geometry.minkowski_reduction.minkowski_reduce`."""
317 from ase.geometry.minkowski_reduction import minkowski_reduce
318 cell, op = minkowski_reduce(self, self.mask())
319 result = Cell(cell)
320 return result, op
322 def permute_axes(self, permutation):
323 """Permute axes of cell."""
324 assert (np.sort(permutation) == np.arange(3)).all()
325 permuted = Cell(self[permutation][:, permutation])
326 return permuted
328 def standard_form(self):
329 """Rotate axes such that unit cell is lower triangular. The cell
330 handedness is preserved.
332 A lower-triangular cell with positive diagonal entries is a canonical
333 (i.e. unique) description. For a left-handed cell the diagonal entries
334 are negative.
336 Returns:
338 rcell: the standardized cell object
340 Q: ndarray
341 The orthogonal transformation. Here, rcell @ Q = cell, where cell
342 is the input cell and rcell is the lower triangular (output) cell.
343 """
345 sign = self.handedness
346 if sign == 0:
347 sign = 1
349 # LQ decomposition provides an axis-aligned description of the cell.
350 # Q is an orthogonal matrix and L is a lower triangular matrix. The
351 # decomposition is a unique description if the diagonal elements are
352 # all positive (negative for a left-handed cell).
353 Q, L = np.linalg.qr(self.T)
354 Q = Q.T
355 L = L.T
357 # correct the signs of the diagonal elements
358 signs = np.sign(np.diag(L))
359 indices = np.where(signs == 0)[0]
360 signs[indices] = 1
361 indices = np.where(signs != sign)[0]
362 L[:, indices] *= -1
363 Q[indices] *= -1
364 return Cell(L), Q
366 # XXX We want a reduction function that brings the cell into
367 # standard form as defined by Setyawan and Curtarolo.