1import numpy as np 

2 

3 

4def cut(atoms, a=(1, 0, 0), b=(0, 1, 0), c=None, clength=None, 

5 origo=(0, 0, 0), nlayers=None, extend=1.0, tolerance=0.01, 

6 maxatoms=None): 

7 """Cuts out a cell defined by *a*, *b*, *c* and *origo* from a 

8 sufficiently repeated copy of *atoms*. 

9 

10 Typically, this function is used to create slabs of different 

11 sizes and orientations. The vectors *a*, *b* and *c* are in scaled 

12 coordinates and defines the returned cell and should normally be 

13 integer-valued in order to end up with a periodic 

14 structure. However, for systems with sub-translations, like fcc, 

15 integer multiples of 1/2 or 1/3 might also make sense for some 

16 directions (and will be treated correctly). 

17 

18 Parameters: 

19 

20 atoms: Atoms instance 

21 This should correspond to a repeatable unit cell. 

22 a: int | 3 floats 

23 The a-vector in scaled coordinates of the cell to cut out. If 

24 integer, the a-vector will be the scaled vector from *origo* to the 

25 atom with index *a*. 

26 b: int | 3 floats 

27 The b-vector in scaled coordinates of the cell to cut out. If 

28 integer, the b-vector will be the scaled vector from *origo* to the 

29 atom with index *b*. 

30 c: None | int | 3 floats 

31 The c-vector in scaled coordinates of the cell to cut out. 

32 if integer, the c-vector will be the scaled vector from *origo* to 

33 the atom with index *c*. 

34 If *None* it will be along cross(a, b) converted to real space 

35 and normalised with the cube root of the volume. Note that this 

36 in general is not perpendicular to a and b for non-cubic 

37 systems. For cubic systems however, this is redused to 

38 c = cross(a, b). 

39 clength: None | float 

40 If not None, the length of the c-vector will be fixed to 

41 *clength* Angstroms. Should not be used together with 

42 *nlayers*. 

43 origo: int | 3 floats 

44 Position of origo of the new cell in scaled coordinates. If 

45 integer, the position of the atom with index *origo* is used. 

46 nlayers: None | int 

47 If *nlayers* is not *None*, the returned cell will have 

48 *nlayers* atomic layers in the c-direction. 

49 extend: 1 or 3 floats 

50 The *extend* argument scales the effective cell in which atoms 

51 will be included. It must either be three floats or a single 

52 float scaling all 3 directions. By setting to a value just 

53 above one, e.g. 1.05, it is possible to all the corner and 

54 edge atoms in the returned cell. This will of cause make the 

55 returned cell non-repeatable, but is very useful for 

56 visualisation. 

57 tolerance: float 

58 Determines what is defined as a plane. All atoms within 

59 *tolerance* Angstroms from a given plane will be considered to 

60 belong to that plane. 

61 maxatoms: None | int 

62 This option is used to auto-tune *tolerance* when *nlayers* is 

63 given for high zone axis systems. For high zone axis one 

64 needs to reduce *tolerance* in order to distinguise the atomic 

65 planes, resulting in the more atoms will be added and 

66 eventually MemoryError. A too small *tolerance*, on the other 

67 hand, might result in inproper splitting of atomic planes and 

68 that too few layers are returned. If *maxatoms* is not None, 

69 *tolerance* will automatically be gradually reduced until 

70 *nlayers* atomic layers is obtained, when the number of atoms 

71 exceeds *maxatoms*. 

72 

73 Example: 

74 

75 >>> import ase 

76 >>> from ase.spacegroup import crystal 

77 >>> 

78 # Create an aluminium (111) slab with three layers 

79 # 

80 # First an unit cell of Al 

81 >>> a = 4.05 

82 >>> aluminium = crystal('Al', [(0,0,0)], spacegroup=225, 

83 ... cellpar=[a, a, a, 90, 90, 90]) 

84 >>> 

85 # Then cut out the slab 

86 >>> al111 = cut(aluminium, (1,-1,0), (0,1,-1), nlayers=3) 

87 >>> 

88 # Visualisation of the skutterudite unit cell 

89 # 

90 # Again, create a skutterudite unit cell 

91 >>> a = 9.04 

92 >>> skutterudite = crystal( 

93 ... ('Co', 'Sb'), 

94 ... basis=[(0.25,0.25,0.25), (0.0, 0.335, 0.158)], 

95 ... spacegroup=204, 

96 ... cellpar=[a, a, a, 90, 90, 90]) 

97 >>> 

98 # Then use *origo* to put 'Co' at the corners and *extend* to 

99 # include all corner and edge atoms. 

100 >>> s = cut(skutterudite, origo=(0.25, 0.25, 0.25), extend=1.01) 

101 >>> ase.view(s) # doctest: +SKIP 

102 """ 

103 atoms = atoms.copy() 

104 cell = atoms.cell 

105 

106 if isinstance(origo, int): 

107 origo = atoms.get_scaled_positions()[origo] 

108 origo = np.array(origo, dtype=float) 

109 

110 scaled = (atoms.get_scaled_positions() - origo) % 1.0 

111 scaled %= 1.0 # needed to ensure that all numbers are *less* than one 

112 atoms.set_scaled_positions(scaled) 

113 

114 if isinstance(a, int): 

115 a = scaled[a] - origo 

116 if isinstance(b, int): 

117 b = scaled[b] - origo 

118 if isinstance(c, int): 

119 c = scaled[c] - origo 

120 

121 a = np.array(a, dtype=float) 

122 b = np.array(b, dtype=float) 

123 if c is None: 

124 metric = np.dot(cell, cell.T) 

125 vol = np.sqrt(np.linalg.det(metric)) 

126 h = np.cross(a, b) 

127 H = np.linalg.solve(metric.T, h.T) 

128 c = vol * H / vol**(1. / 3.) 

129 c = np.array(c, dtype=float) 

130 

131 if nlayers: 

132 # Recursive increase the length of c until we have at least 

133 # *nlayers* atomic layers parallel to the a-b plane 

134 while True: 

135 at = cut(atoms, a, b, c, origo=origo, extend=extend, 

136 tolerance=tolerance) 

137 scaled = at.get_scaled_positions() 

138 d = scaled[:, 2] 

139 keys = np.argsort(d) 

140 ikeys = np.argsort(keys) 

141 tol = tolerance 

142 while True: 

143 mask = np.concatenate(([True], np.diff(d[keys]) > tol)) 

144 tags = np.cumsum(mask)[ikeys] - 1 

145 levels = d[keys][mask] 

146 if (maxatoms is None or len(at) < maxatoms or 

147 len(levels) > nlayers): 

148 break 

149 tol *= 0.9 

150 if len(levels) > nlayers: 

151 break 

152 c *= 2 

153 

154 at.cell[2] *= levels[nlayers] 

155 return at[tags < nlayers] 

156 

157 newcell = np.dot(np.array([a, b, c]), cell) 

158 if nlayers is None and clength is not None: 

159 newcell[2, :] *= clength / np.linalg.norm(newcell[2]) 

160 

161 # Create a new atoms object, repeated and translated such that 

162 # it completely covers the new cell 

163 scorners_newcell = np.array([[0., 0., 0.], [0., 0., 1.], 

164 [0., 1., 0.], [0., 1., 1.], 

165 [1., 0., 0.], [1., 0., 1.], 

166 [1., 1., 0.], [1., 1., 1.]]) 

167 corners = np.dot(scorners_newcell, newcell * extend) 

168 scorners = np.linalg.solve(cell.T, corners.T).T 

169 rep = np.ceil(scorners.ptp(axis=0)).astype('int') + 1 

170 trans = np.dot(np.floor(scorners.min(axis=0)), cell) 

171 atoms = atoms.repeat(rep) 

172 atoms.translate(trans) 

173 atoms.set_cell(newcell) 

174 

175 # Mask out atoms outside new cell 

176 stol = 0.1 * tolerance # scaled tolerance, XXX 

177 maskcell = atoms.cell * extend 

178 sp = np.linalg.solve(maskcell.T, (atoms.positions).T).T 

179 mask = np.all(np.logical_and(-stol <= sp, sp < 1 - stol), axis=1) 

180 atoms = atoms[mask] 

181 return atoms 

182 

183 

184class IncompatibleCellError(ValueError): 

185 """Exception raised if stacking fails due to incompatible cells 

186 between *atoms1* and *atoms2*.""" 

187 pass 

188 

189 

190def stack(atoms1, atoms2, axis=2, cell=None, fix=0.5, 

191 maxstrain=0.5, distance=None, reorder=False, 

192 output_strained=False): 

193 """Return a new Atoms instance with *atoms2* stacked on top of 

194 *atoms1* along the given axis. Periodicity in all directions is 

195 ensured. 

196 

197 The size of the final cell is determined by *cell*, except 

198 that the length alongh *axis* will be the sum of 

199 *atoms1.cell[axis]* and *atoms2.cell[axis]*. If *cell* is None, 

200 it will be interpolated between *atoms1* and *atoms2*, where 

201 *fix* determines their relative weight. Hence, if *fix* equals 

202 zero, the final cell will be determined purely from *atoms1* and 

203 if *fix* equals one, it will be determined purely from 

204 *atoms2*. 

205 

206 An ase.geometry.IncompatibleCellError exception is raised if the 

207 cells of *atoms1* and *atoms2* are incompatible, e.g. if the far 

208 corner of the unit cell of either *atoms1* or *atoms2* is 

209 displaced more than *maxstrain*. Setting *maxstrain* to None 

210 disables this check. 

211 

212 If *distance* is not None, the size of the final cell, along the 

213 direction perpendicular to the interface, will be adjusted such 

214 that the distance between the closest atoms in *atoms1* and 

215 *atoms2* will be equal to *distance*. This option uses 

216 scipy.optimize.fmin() and hence require scipy to be installed. 

217 

218 If *reorder* is True, then the atoms will be reordered such that 

219 all atoms with the same symbol will follow sequencially after each 

220 other, eg: 'Al2MnAl10Fe' -> 'Al12FeMn'. 

221 

222 If *output_strained* is True, then the strained versions of 

223 *atoms1* and *atoms2* are returned in addition to the stacked 

224 structure. 

225 

226 Example: 

227 

228 >>> import ase 

229 >>> from ase.spacegroup import crystal 

230 >>> 

231 # Create an Ag(110)-Si(110) interface with three atomic layers 

232 # on each side. 

233 >>> a_ag = 4.09 

234 >>> ag = crystal(['Ag'], basis=[(0,0,0)], spacegroup=225, 

235 ... cellpar=[a_ag, a_ag, a_ag, 90., 90., 90.]) 

236 >>> ag110 = cut(ag, (0, 0, 3), (-1.5, 1.5, 0), nlayers=3) 

237 >>> 

238 >>> a_si = 5.43 

239 >>> si = crystal(['Si'], basis=[(0,0,0)], spacegroup=227, 

240 ... cellpar=[a_si, a_si, a_si, 90., 90., 90.]) 

241 >>> si110 = cut(si, (0, 0, 2), (-1, 1, 0), nlayers=3) 

242 >>> 

243 >>> interface = stack(ag110, si110, maxstrain=1) 

244 >>> ase.view(interface) # doctest: +SKIP 

245 >>> 

246 # Once more, this time adjusted such that the distance between 

247 # the closest Ag and Si atoms will be 2.3 Angstrom (requires scipy). 

248 >>> interface2 = stack(ag110, si110, 

249 ... maxstrain=1, distance=2.3) # doctest:+ELLIPSIS 

250 Optimization terminated successfully. 

251 ... 

252 >>> ase.view(interface2) # doctest: +SKIP 

253 """ 

254 atoms1 = atoms1.copy() 

255 atoms2 = atoms2.copy() 

256 

257 for atoms in [atoms1, atoms2]: 

258 if not atoms.cell[axis].any(): 

259 atoms.center(vacuum=0.0, axis=axis) 

260 

261 if (np.sign(np.linalg.det(atoms1.cell)) != 

262 np.sign(np.linalg.det(atoms2.cell))): 

263 raise IncompatibleCellError('Cells of *atoms1* and *atoms2* must have ' 

264 'same handedness.') 

265 

266 c1 = np.linalg.norm(atoms1.cell[axis]) 

267 c2 = np.linalg.norm(atoms2.cell[axis]) 

268 if cell is None: 

269 cell1 = atoms1.cell.copy() 

270 cell2 = atoms2.cell.copy() 

271 cell1[axis] /= c1 

272 cell2[axis] /= c2 

273 cell = cell1 + fix * (cell2 - cell1) 

274 cell[axis] /= np.linalg.norm(cell[axis]) 

275 cell1 = cell.copy() 

276 cell2 = cell.copy() 

277 cell1[axis] *= c1 

278 cell2[axis] *= c2 

279 

280 if maxstrain: 

281 strain1 = np.sqrt(((cell1 - atoms1.cell).sum(axis=0)**2).sum()) 

282 strain2 = np.sqrt(((cell2 - atoms2.cell).sum(axis=0)**2).sum()) 

283 if strain1 > maxstrain or strain2 > maxstrain: 

284 raise IncompatibleCellError( 

285 '*maxstrain* exceeded. *atoms1* strained %f and ' 

286 '*atoms2* strained %f.' % (strain1, strain2)) 

287 

288 atoms1.set_cell(cell1, scale_atoms=True) 

289 atoms2.set_cell(cell2, scale_atoms=True) 

290 if output_strained: 

291 atoms1_strained = atoms1.copy() 

292 atoms2_strained = atoms2.copy() 

293 

294 if distance is not None: 

295 from scipy.optimize import fmin 

296 

297 def mindist(pos1, pos2): 

298 n1 = len(pos1) 

299 n2 = len(pos2) 

300 idx1 = np.arange(n1).repeat(n2) 

301 idx2 = np.tile(np.arange(n2), n1) 

302 return np.sqrt(((pos1[idx1] - pos2[idx2])**2).sum(axis=1).min()) 

303 

304 def func(x): 

305 t1, t2, h1, h2 = x[0:3], x[3:6], x[6], x[7] 

306 pos1 = atoms1.positions + t1 

307 pos2 = atoms2.positions + t2 

308 d1 = mindist(pos1, pos2 + (h1 + 1.0) * atoms1.cell[axis]) 

309 d2 = mindist(pos2, pos1 + (h2 + 1.0) * atoms2.cell[axis]) 

310 return (d1 - distance)**2 + (d2 - distance)**2 

311 

312 atoms1.center() 

313 atoms2.center() 

314 x0 = np.zeros((8,)) 

315 x = fmin(func, x0) 

316 t1, t2, h1, h2 = x[0:3], x[3:6], x[6], x[7] 

317 atoms1.translate(t1) 

318 atoms2.translate(t2) 

319 atoms1.cell[axis] *= 1.0 + h1 

320 atoms2.cell[axis] *= 1.0 + h2 

321 

322 atoms2.translate(atoms1.cell[axis]) 

323 atoms1.cell[axis] += atoms2.cell[axis] 

324 atoms1.extend(atoms2) 

325 

326 if reorder: 

327 atoms1 = sort(atoms1) 

328 

329 if output_strained: 

330 return atoms1, atoms1_strained, atoms2_strained 

331 else: 

332 return atoms1 

333 

334 

335def rotation_matrix(a1, a2, b1, b2): 

336 """Returns a rotation matrix that rotates the vectors *a1* in the 

337 direction of *a2* and *b1* in the direction of *b2*. 

338 

339 In the case that the angle between *a2* and *b2* is not the same 

340 as between *a1* and *b1*, a proper rotation matrix will anyway be 

341 constructed by first rotate *b2* in the *b1*, *b2* plane. 

342 """ 

343 a1 = np.asarray(a1, dtype=float) / np.linalg.norm(a1) 

344 b1 = np.asarray(b1, dtype=float) / np.linalg.norm(b1) 

345 c1 = np.cross(a1, b1) 

346 c1 /= np.linalg.norm(c1) # clean out rounding errors... 

347 

348 a2 = np.asarray(a2, dtype=float) / np.linalg.norm(a2) 

349 b2 = np.asarray(b2, dtype=float) / np.linalg.norm(b2) 

350 c2 = np.cross(a2, b2) 

351 c2 /= np.linalg.norm(c2) # clean out rounding errors... 

352 

353 # Calculate rotated *b2* 

354 theta = np.arccos(np.dot(a2, b2)) - np.arccos(np.dot(a1, b1)) 

355 b3 = np.sin(theta) * a2 + np.cos(theta) * b2 

356 b3 /= np.linalg.norm(b3) # clean out rounding errors... 

357 

358 A1 = np.array([a1, b1, c1]) 

359 A2 = np.array([a2, b3, c2]) 

360 R = np.linalg.solve(A1, A2).T 

361 return R 

362 

363 

364def rotate(atoms, a1, a2, b1, b2, rotate_cell=True, center=(0, 0, 0)): 

365 """Rotate *atoms*, such that *a1* will be rotated in the direction 

366 of *a2* and *b1* in the direction of *b2*. The point at *center* 

367 is fixed. Use *center='COM'* to fix the center of mass. If 

368 *rotate_cell* is true, the cell will be rotated together with the 

369 atoms. 

370 

371 Note that the 000-corner of the cell is by definition fixed at 

372 origo. Hence, setting *center* to something other than (0, 0, 0) 

373 will rotate the atoms out of the cell, even if *rotate_cell* is 

374 True. 

375 """ 

376 if isinstance(center, str) and center.lower() == 'com': 

377 center = atoms.get_center_of_mass() 

378 

379 R = rotation_matrix(a1, a2, b1, b2) 

380 atoms.positions[:] = np.dot(atoms.positions - center, R.T) + center 

381 

382 if rotate_cell: 

383 atoms.cell[:] = np.dot(atoms.cell, R.T) 

384 

385 

386def minimize_tilt_ij(atoms, modified=1, fixed=0, fold_atoms=True): 

387 """Minimize the tilt angle for two given axes. 

388 

389 The problem is underdetermined. Therefore one can choose one axis 

390 that is kept fixed. 

391 """ 

392 

393 orgcell_cc = atoms.get_cell() 

394 pbc_c = atoms.get_pbc() 

395 i = fixed 

396 j = modified 

397 if not (pbc_c[i] and pbc_c[j]): 

398 raise RuntimeError('Axes have to be periodic') 

399 

400 prod_cc = np.dot(orgcell_cc, orgcell_cc.T) 

401 cell_cc = 1. * orgcell_cc 

402 nji = np.floor(- prod_cc[i, j] / prod_cc[i, i] + 0.5) 

403 cell_cc[j] = orgcell_cc[j] + nji * cell_cc[i] 

404 

405 # sanity check 

406 def volume(cell): 

407 return np.abs(np.dot(cell[2], np.cross(cell[0], cell[1]))) 

408 V = volume(cell_cc) 

409 assert(abs(volume(orgcell_cc) - V) / V < 1.e-10) 

410 

411 atoms.set_cell(cell_cc) 

412 

413 if fold_atoms: 

414 atoms.wrap() 

415 

416 

417def minimize_tilt(atoms, order=range(3), fold_atoms=True): 

418 """Minimize the tilt angles of the unit cell.""" 

419 pbc_c = atoms.get_pbc() 

420 

421 for i1, c1 in enumerate(order): 

422 for c2 in order[i1 + 1:]: 

423 if pbc_c[c1] and pbc_c[c2]: 

424 minimize_tilt_ij(atoms, c1, c2, fold_atoms) 

425 

426 

427def niggli_reduce_cell(cell, epsfactor=None): 

428 from ase.geometry import cellpar_to_cell 

429 

430 if epsfactor is None: 

431 epsfactor = 1e-5 

432 eps = epsfactor * abs(np.linalg.det(cell))**(1./3.) 

433 

434 cell = np.asarray(cell) 

435 

436 I3 = np.eye(3, dtype=int) 

437 I6 = np.eye(6, dtype=int) 

438 

439 C = I3.copy() 

440 D = I6.copy() 

441 

442 g0 = np.zeros(6, dtype=float) 

443 g0[0] = np.dot(cell[0], cell[0]) 

444 g0[1] = np.dot(cell[1], cell[1]) 

445 g0[2] = np.dot(cell[2], cell[2]) 

446 g0[3] = 2 * np.dot(cell[1], cell[2]) 

447 g0[4] = 2 * np.dot(cell[0], cell[2]) 

448 g0[5] = 2 * np.dot(cell[0], cell[1]) 

449 

450 g = np.dot(D, g0) 

451 

452 def lt(x, y, eps=eps): 

453 return x < y - eps 

454 

455 def gt(x, y, eps=eps): 

456 return lt(y, x, eps) 

457 

458 def eq(x, y, eps=eps): 

459 return not (lt(x, y, eps) or gt(x, y, eps)) 

460 

461 for _ in range(10000): 

462 if (gt(g[0], g[1]) 

463 or (eq(g[0], g[1]) and gt(abs(g[3]), abs(g[4])))): 

464 C = np.dot(C, -I3[[1, 0, 2]]) 

465 D = np.dot(I6[[1, 0, 2, 4, 3, 5]], D) 

466 g = np.dot(D, g0) 

467 continue 

468 elif (gt(g[1], g[2]) 

469 or (eq(g[1], g[2]) and gt(abs(g[4]), abs(g[5])))): 

470 C = np.dot(C, -I3[[0, 2, 1]]) 

471 D = np.dot(I6[[0, 2, 1, 3, 5, 4]], D) 

472 g = np.dot(D, g0) 

473 continue 

474 

475 lmn = np.array(gt(g[3:], 0, eps=eps/2), dtype=int) 

476 lmn -= np.array(lt(g[3:], 0, eps=eps/2), dtype=int) 

477 

478 if lmn.prod() == 1: 

479 ijk = lmn.copy() 

480 for idx in range(3): 

481 if ijk[idx] == 0: 

482 ijk[idx] = 1 

483 else: 

484 ijk = np.ones(3, dtype=int) 

485 if np.any(lmn != -1): 

486 r = None 

487 for idx in range(3): 

488 if lmn[idx] == 1: 

489 ijk[idx] = -1 

490 elif lmn[idx] == 0: 

491 r = idx 

492 if ijk.prod() == -1: 

493 ijk[r] = -1 

494 

495 C *= ijk[np.newaxis] 

496 

497 D[3] *= ijk[1] * ijk[2] 

498 D[4] *= ijk[0] * ijk[2] 

499 D[5] *= ijk[0] * ijk[1] 

500 g = np.dot(D, g0) 

501 

502 if (gt(abs(g[3]), g[1]) 

503 or (eq(g[3], g[1]) and lt(2 * g[4], g[5])) 

504 or (eq(g[3], -g[1]) and lt(g[5], 0))): 

505 s = int(np.sign(g[3])) 

506 

507 A = I3.copy() 

508 A[1, 2] = -s 

509 C = np.dot(C, A) 

510 

511 B = I6.copy() 

512 B[2, 1] = 1 

513 B[2, 3] = -s 

514 B[3, 1] = -2 * s 

515 B[4, 5] = -s 

516 D = np.dot(B, D) 

517 g = np.dot(D, g0) 

518 elif (gt(abs(g[4]), g[0]) 

519 or (eq(g[4], g[0]) and lt(2 * g[3], g[5])) 

520 or (eq(g[4], -g[0]) and lt(g[5], 0))): 

521 s = int(np.sign(g[4])) 

522 

523 A = I3.copy() 

524 A[0, 2] = -s 

525 C = np.dot(C, A) 

526 

527 B = I6.copy() 

528 B[2, 0] = 1 

529 B[2, 4] = -s 

530 B[3, 5] = -s 

531 B[4, 0] = -2 * s 

532 D = np.dot(B, D) 

533 g = np.dot(D, g0) 

534 elif (gt(abs(g[5]), g[0]) 

535 or (eq(g[5], g[0]) and lt(2 * g[3], g[4])) 

536 or (eq(g[5], -g[0]) and lt(g[4], 0))): 

537 s = int(np.sign(g[5])) 

538 

539 A = I3.copy() 

540 A[0, 1] = -s 

541 C = np.dot(C, A) 

542 

543 B = I6.copy() 

544 B[1, 0] = 1 

545 B[1, 5] = -s 

546 B[3, 4] = -s 

547 B[5, 0] = -2 * s 

548 D = np.dot(B, D) 

549 g = np.dot(D, g0) 

550 elif (lt(g[[0, 1, 3, 4, 5]].sum(), 0) 

551 or (eq(g[[0, 1, 3, 4, 5]].sum(), 0) 

552 and gt(2 * (g[0] + g[4]) + g[5], 0))): 

553 A = I3.copy() 

554 A[:, 2] = 1 

555 C = np.dot(C, A) 

556 

557 B = I6.copy() 

558 B[2, :] = 1 

559 B[3, 1] = 2 

560 B[3, 5] = 1 

561 B[4, 0] = 2 

562 B[4, 5] = 1 

563 D = np.dot(B, D) 

564 g = np.dot(D, g0) 

565 else: 

566 break 

567 else: 

568 raise RuntimeError('Niggli reduction not done in 10000 steps!\n' 

569 'cell={}\n' 

570 'operation={}' 

571 .format(cell.tolist(), C.tolist())) 

572 

573 abc = np.sqrt(g[:3]) 

574 # Prevent division by zero e.g. for cell==zeros((3, 3)): 

575 abcprod = max(abc.prod(), 1e-100) 

576 cosangles = abc * g[3:] / (2 * abcprod) 

577 angles = 180 * np.arccos(cosangles) / np.pi 

578 newcell = np.array(cellpar_to_cell(np.concatenate([abc, angles])), 

579 dtype=float) 

580 

581 return newcell, C 

582 

583 

584def update_cell_and_positions(atoms, new_cell, op): 

585 """Helper method for transforming cell and positions of atoms object.""" 

586 scpos = np.linalg.solve(op, atoms.get_scaled_positions().T).T 

587 scpos %= 1.0 

588 scpos %= 1.0 

589 

590 atoms.set_cell(new_cell) 

591 atoms.set_scaled_positions(scpos) 

592 

593 

594def niggli_reduce(atoms): 

595 """Convert the supplied atoms object's unit cell into its 

596 maximally-reduced Niggli unit cell. Even if the unit cell is already 

597 maximally reduced, it will be converted into its unique Niggli unit cell. 

598 This will also wrap all atoms into the new unit cell. 

599 

600 References: 

601 

602 Niggli, P. "Krystallographische und strukturtheoretische Grundbegriffe. 

603 Handbuch der Experimentalphysik", 1928, Vol. 7, Part 1, 108-176. 

604 

605 Krivy, I. and Gruber, B., "A Unified Algorithm for Determining the 

606 Reduced (Niggli) Cell", Acta Cryst. 1976, A32, 297-298. 

607 

608 Grosse-Kunstleve, R.W.; Sauter, N. K.; and Adams, P. D. "Numerically 

609 stable algorithms for the computation of reduced unit cells", Acta Cryst. 

610 2004, A60, 1-6. 

611 """ 

612 

613 assert all(atoms.pbc), 'Can only reduce 3d periodic unit cells!' 

614 new_cell, op = niggli_reduce_cell(atoms.cell) 

615 update_cell_and_positions(atoms, new_cell, op) 

616 

617 

618def reduce_lattice(atoms, eps=2e-4): 

619 """Reduce atoms object to canonical lattice. 

620 

621 This changes the cell and positions such that the atoms object has 

622 the canonical form used for defining band paths but is otherwise 

623 physically equivalent. The eps parameter is used as a tolerance 

624 for determining the cell's Bravais lattice.""" 

625 from ase.geometry.bravais_type_engine import identify_lattice 

626 niggli_reduce(atoms) 

627 lat, op = identify_lattice(atoms.cell, eps=eps) 

628 update_cell_and_positions(atoms, lat.tocell(), np.linalg.inv(op)) 

629 

630 

631def sort(atoms, tags=None): 

632 """Return a new Atoms object with sorted atomic order. The default 

633 is to order according to chemical symbols, but if *tags* is not 

634 None, it will be used instead. A stable sorting algorithm is used. 

635 

636 Example: 

637 

638 >>> from ase.build import bulk 

639 >>> # Two unit cells of NaCl: 

640 >>> a = 5.64 

641 >>> nacl = bulk('NaCl', 'rocksalt', a=a) * (2, 1, 1) 

642 >>> nacl.get_chemical_symbols() 

643 ['Na', 'Cl', 'Na', 'Cl'] 

644 >>> nacl_sorted = sort(nacl) 

645 >>> nacl_sorted.get_chemical_symbols() 

646 ['Cl', 'Cl', 'Na', 'Na'] 

647 >>> np.all(nacl_sorted.cell == nacl.cell) 

648 True 

649 """ 

650 if tags is None: 

651 tags = atoms.get_chemical_symbols() 

652 else: 

653 tags = list(tags) 

654 deco = sorted([(tag, i) for i, tag in enumerate(tags)]) 

655 indices = [i for tag, i in deco] 

656 return atoms[indices]