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.