1import ase 

2from typing import Mapping, Sequence, Union 

3import numpy as np 

4from ase.utils.arraywrapper import arraylike 

5from ase.utils import pbc2pbc 

6 

7 

8__all__ = ['Cell'] 

9 

10 

11@arraylike 

12class Cell: 

13 """Parallel epipedal unit cell of up to three dimensions. 

14 

15 This object resembles a 3x3 array whose [i, j]-th element is the jth 

16 Cartesian coordinate of the ith unit vector. 

17 

18 Cells of less than three dimensions are represented by placeholder 

19 unit vectors that are zero.""" 

20 

21 ase_objtype = 'cell' # For JSON'ing 

22 

23 def __init__(self, array): 

24 """Create cell. 

25 

26 Parameters: 

27 

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 

34 

35 def cellpar(self, radians=False): 

36 """Get unit cell parameters. Sequence of 6 numbers. 

37 

38 First three are unit cell vector lengths and second three 

39 are angles between them:: 

40 

41 [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)] 

42 

43 in degrees. 

44 

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) 

48 

49 def todict(self): 

50 return dict(array=self.array) 

51 

52 @classmethod 

53 def ascell(cls, cell): 

54 """Return argument as a Cell object. See :meth:`ase.cell.Cell.new`. 

55 

56 A new Cell object is created if necessary.""" 

57 if isinstance(cell, cls): 

58 return cell 

59 return cls.new(cell) 

60 

61 @classmethod 

62 def new(cls, cell=None): 

63 """Create new cell from any parameters. 

64 

65 If cell is three numbers, assume three lengths with right angles. 

66 

67 If cell is six numbers, assume three lengths, then three angles. 

68 

69 If cell is 3x3, assume three cell vectors.""" 

70 

71 if cell is None: 

72 cell = np.zeros((3, 3)) 

73 

74 cell = np.array(cell, float) 

75 

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!') 

84 

85 cellobj = cls(cell) 

86 return cellobj 

87 

88 @classmethod 

89 def fromcellpar(cls, cellpar, ab_normal=(0, 0, 1), a_direction=None): 

90 """Return new Cell from cell lengths and angles. 

91 

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) 

96 

97 def get_bravais_lattice(self, eps=2e-4, *, pbc=True): 

98 """Return :class:`~ase.lattice.BravaisLattice` for this cell: 

99 

100 >>> cell = Cell.fromcellpar([4, 4, 4, 60, 60, 60]) 

101 >>> print(cell.get_bravais_lattice()) 

102 FCC(a=5.65685) 

103 

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. 

109 

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. 

113 

114 """ 

115 from ase.lattice import identify_lattice 

116 pbc = self.any(1) & pbc2pbc(pbc) 

117 lat, op = identify_lattice(self, eps=eps, pbc=pbc) 

118 return lat 

119 

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. 

131 

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. 

135 

136 If special special points are given, interpolate the path 

137 directly from the available data. 

138 

139 Parameters: 

140 

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. 

157 

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]) 

163 

164 """ 

165 # TODO: Combine with the rotation transformation from bandpath() 

166 

167 cell = self.uncomplete(pbc) 

168 

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) 

180 

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 

188 

189 def complete(self): 

190 """Convert missing cell vectors into orthogonal unit vectors.""" 

191 from ase.geometry.cell import complete_cell 

192 cell = Cell(complete_cell(self.array)) 

193 return cell 

194 

195 def copy(self): 

196 """Return a copy of this cell.""" 

197 cell = Cell(self.array.copy()) 

198 return cell 

199 

200 @property 

201 def rank(self) -> int: 

202 """"Return the dimension of the cell. 

203 

204 Equal to the number of nonzero lattice vectors.""" 

205 # The name ndim clashes with ndarray.ndim 

206 return self.any(1).sum() # type: ignore 

207 

208 @property 

209 def orthorhombic(self) -> bool: 

210 """Return whether this cell is represented by a diagonal matrix.""" 

211 from ase.geometry.cell import is_orthorhombic 

212 return is_orthorhombic(self) 

213 

214 def lengths(self): 

215 """Return the length of each lattice vector as an array.""" 

216 return np.linalg.norm(self, axis=1) 

217 

218 def angles(self): 

219 """Return an array with the three angles alpha, beta, and gamma.""" 

220 return self.cellpar()[3:].copy() 

221 

222 def __array__(self, dtype=float): 

223 if dtype != float: 

224 raise ValueError('Cannot convert cell to array of type {}' 

225 .format(dtype)) 

226 return self.array 

227 

228 def __bool__(self): 

229 return bool(self.any()) # need to convert from np.bool_ 

230 

231 __nonzero__ = __bool__ 

232 

233 @property 

234 def volume(self) -> float: 

235 """Get the volume of this cell. 

236 

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)) 

242 

243 @property 

244 def handedness(self) -> int: 

245 """Sign of the determinant of the matrix of cell vectors. 

246 

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))) 

250 

251 def scaled_positions(self, positions) -> np.ndarray: 

252 """Calculate scaled positions from Cartesian positions. 

253 

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 

259 

260 def cartesian_positions(self, scaled_positions) -> np.ndarray: 

261 """Calculate Cartesian positions from scaled positions.""" 

262 return scaled_positions @ self.complete() 

263 

264 def reciprocal(self) -> 'Cell': 

265 """Get reciprocal lattice as a Cell object. 

266 

267 Does not include factor of 2 pi.""" 

268 return Cell(np.linalg.pinv(self).transpose()) 

269 

270 def __repr__(self): 

271 if self.orthorhombic: 

272 numbers = self.lengths().tolist() 

273 else: 

274 numbers = self.tolist() 

275 

276 return 'Cell({})'.format(numbers) 

277 

278 def niggli_reduce(self, eps=1e-5): 

279 """Niggli reduce this cell, returning a new cell and mapping. 

280 

281 See also :func:`ase.build.tools.niggli_reduce_cell`.""" 

282 from ase.build.tools import niggli_reduce_cell 

283 cell, op = niggli_reduce_cell(self, epsfactor=eps) 

284 result = Cell(cell) 

285 return result, op 

286 

287 def minkowski_reduce(self): 

288 """Minkowski-reduce this cell, returning new cell and mapping. 

289 

290 See also :func:`ase.geometry.minkowski_reduction.minkowski_reduce`.""" 

291 from ase.geometry.minkowski_reduction import minkowski_reduce 

292 cell, op = minkowski_reduce(self, self.any(1)) 

293 result = Cell(cell) 

294 return result, op 

295 

296 def permute_axes(self, permutation): 

297 """Permute axes of cell.""" 

298 assert (np.sort(permutation) == np.arange(3)).all() 

299 permuted = Cell(self[permutation][:, permutation]) 

300 return permuted 

301 

302 def standard_form(self): 

303 """Rotate axes such that unit cell is lower triangular. The cell 

304 handedness is preserved. 

305 

306 A lower-triangular cell with positive diagonal entries is a canonical 

307 (i.e. unique) description. For a left-handed cell the diagonal entries 

308 are negative. 

309 

310 Returns: 

311 

312 rcell: the standardized cell object 

313 

314 Q: ndarray 

315 The orthogonal transformation. Here, rcell @ Q = cell, where cell 

316 is the input cell and rcell is the lower triangular (output) cell. 

317 """ 

318 

319 # get cell handedness (right or left) 

320 sign = np.sign(np.linalg.det(self)) 

321 if sign == 0: 

322 sign = 1 

323 

324 # LQ decomposition provides an axis-aligned description of the cell. 

325 # Q is an orthogonal matrix and L is a lower triangular matrix. The 

326 # decomposition is a unique description if the diagonal elements are 

327 # all positive (negative for a left-handed cell). 

328 Q, L = np.linalg.qr(self.T) 

329 Q = Q.T 

330 L = L.T 

331 

332 # correct the signs of the diagonal elements 

333 signs = np.sign(np.diag(L)) 

334 indices = np.where(signs == 0)[0] 

335 signs[indices] = 1 

336 indices = np.where(signs != sign)[0] 

337 L[:, indices] *= -1 

338 Q[indices] *= -1 

339 return Cell(L), Q 

340 

341 # XXX We want a reduction function that brings the cell into 

342 # standard form as defined by Setyawan and Curtarolo.