1import numpy as np 

2from ase.build.niggli import niggli_reduce_cell 

3 

4 

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

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

7 maxatoms=None): 

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

9 sufficiently repeated copy of *atoms*. 

10 

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

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

13 coordinates and defines the returned cell and should normally be 

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

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

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

17 directions (and will be treated correctly). 

18 

19 Parameters: 

20 

21 atoms: Atoms instance 

22 This should correspond to a repeatable unit cell. 

23 a: int | 3 floats 

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

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

26 atom with index *a*. 

27 b: int | 3 floats 

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

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

30 atom with index *b*. 

31 c: None | int | 3 floats 

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

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

34 the atom with index *c*. 

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

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

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

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

39 c = cross(a, b). 

40 clength: None | float 

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

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

43 *nlayers*. 

44 origo: int | 3 floats 

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

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

47 nlayers: None | int 

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

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

50 extend: 1 or 3 floats 

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

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

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

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

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

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

57 visualisation. 

58 tolerance: float 

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

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

61 belong to that plane. 

62 maxatoms: None | int 

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

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

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

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

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

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

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

70 *tolerance* will automatically be gradually reduced until 

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

72 exceeds *maxatoms*. 

73 

74 Example: 

75 

76 >>> import ase 

77 >>> from ase.spacegroup import crystal 

78 >>> 

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

80 # 

81 # First an unit cell of Al 

82 >>> a = 4.05 

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

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

85 >>> 

86 # Then cut out the slab 

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

88 >>> 

89 # Visualisation of the skutterudite unit cell 

90 # 

91 # Again, create a skutterudite unit cell 

92 >>> a = 9.04 

93 >>> skutterudite = crystal( 

94 ... ('Co', 'Sb'), 

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

96 ... spacegroup=204, 

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

98 >>> 

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

100 # include all corner and edge atoms. 

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

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

103 """ 

104 atoms = atoms.copy() 

105 cell = atoms.cell 

106 

107 if isinstance(origo, int): 

108 origo = atoms.get_scaled_positions()[origo] 

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

110 

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

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

113 atoms.set_scaled_positions(scaled) 

114 

115 if isinstance(a, int): 

116 a = scaled[a] - origo 

117 if isinstance(b, int): 

118 b = scaled[b] - origo 

119 if isinstance(c, int): 

120 c = scaled[c] - origo 

121 

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

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

124 if c is None: 

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

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

127 h = np.cross(a, b) 

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

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

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

131 

132 if nlayers: 

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

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

135 while True: 

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

137 tolerance=tolerance) 

138 scaled = at.get_scaled_positions() 

139 d = scaled[:, 2] 

140 keys = np.argsort(d) 

141 ikeys = np.argsort(keys) 

142 tol = tolerance 

143 while True: 

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

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

146 levels = d[keys][mask] 

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

148 len(levels) > nlayers): 

149 break 

150 tol *= 0.9 

151 if len(levels) > nlayers: 

152 break 

153 c *= 2 

154 

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

156 return at[tags < nlayers] 

157 

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

159 if nlayers is None and clength is not None: 

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

161 

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

163 # it completely covers the new cell 

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

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

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

167 [1., 1., 0.], [1., 1., 1.]]) 

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

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

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

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

172 atoms = atoms.repeat(rep) 

173 atoms.translate(trans) 

174 atoms.set_cell(newcell) 

175 

176 # Mask out atoms outside new cell 

177 stol = 0.1 * tolerance # scaled tolerance, XXX 

178 maskcell = atoms.cell * extend 

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

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

181 atoms = atoms[mask] 

182 return atoms 

183 

184 

185class IncompatibleCellError(ValueError): 

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

187 between *atoms1* and *atoms2*.""" 

188 pass 

189 

190 

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

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

193 output_strained=False): 

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

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

196 ensured. 

197 

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

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

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

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

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

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

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

205 *atoms2*. 

206 

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

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

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

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

211 disables this check. 

212 

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

214 direction perpendicular to the interface, will be adjusted such 

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

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

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

218 

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

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

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

222 

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

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

225 structure. 

226 

227 Example: 

228 

229 >>> import ase 

230 >>> from ase.spacegroup import crystal 

231 >>> 

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

233 # on each side. 

234 >>> a_ag = 4.09 

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

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

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

238 >>> 

239 >>> a_si = 5.43 

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

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

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

243 >>> 

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

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

246 >>> 

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

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

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

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

251 Optimization terminated successfully. 

252 ... 

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

254 """ 

255 atoms1 = atoms1.copy() 

256 atoms2 = atoms2.copy() 

257 

258 for atoms in [atoms1, atoms2]: 

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

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

261 

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

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

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

265 'same handedness.') 

266 

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

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

269 if cell is None: 

270 cell1 = atoms1.cell.copy() 

271 cell2 = atoms2.cell.copy() 

272 cell1[axis] /= c1 

273 cell2[axis] /= c2 

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

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

276 cell1 = cell.copy() 

277 cell2 = cell.copy() 

278 cell1[axis] *= c1 

279 cell2[axis] *= c2 

280 

281 if maxstrain: 

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

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

284 if strain1 > maxstrain or strain2 > maxstrain: 

285 raise IncompatibleCellError( 

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

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

288 

289 atoms1.set_cell(cell1, scale_atoms=True) 

290 atoms2.set_cell(cell2, scale_atoms=True) 

291 if output_strained: 

292 atoms1_strained = atoms1.copy() 

293 atoms2_strained = atoms2.copy() 

294 

295 if distance is not None: 

296 from scipy.optimize import fmin 

297 

298 def mindist(pos1, pos2): 

299 n1 = len(pos1) 

300 n2 = len(pos2) 

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

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

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

304 

305 def func(x): 

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

307 pos1 = atoms1.positions + t1 

308 pos2 = atoms2.positions + t2 

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

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

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

312 

313 atoms1.center() 

314 atoms2.center() 

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

316 x = fmin(func, x0) 

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

318 atoms1.translate(t1) 

319 atoms2.translate(t2) 

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

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

322 

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

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

325 atoms1.extend(atoms2) 

326 

327 if reorder: 

328 atoms1 = sort(atoms1) 

329 

330 if output_strained: 

331 return atoms1, atoms1_strained, atoms2_strained 

332 else: 

333 return atoms1 

334 

335 

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

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

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

339 

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

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

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

343 """ 

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

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

346 c1 = np.cross(a1, b1) 

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

348 

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

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

351 c2 = np.cross(a2, b2) 

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

353 

354 # Calculate rotated *b2* 

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

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

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

358 

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

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

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

362 return R 

363 

364 

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

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

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

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

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

370 atoms. 

371 

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

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

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

375 True. 

376 """ 

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

378 center = atoms.get_center_of_mass() 

379 

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

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

382 

383 if rotate_cell: 

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

385 

386 

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

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

389 

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

391 that is kept fixed. 

392 """ 

393 

394 orgcell_cc = atoms.get_cell() 

395 pbc_c = atoms.get_pbc() 

396 i = fixed 

397 j = modified 

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

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

400 

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

402 cell_cc = 1. * orgcell_cc 

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

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

405 

406 # sanity check 

407 def volume(cell): 

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

409 V = volume(cell_cc) 

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

411 

412 atoms.set_cell(cell_cc) 

413 

414 if fold_atoms: 

415 atoms.wrap() 

416 

417 

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

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

420 pbc_c = atoms.get_pbc() 

421 

422 for i1, c1 in enumerate(order): 

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

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

425 minimize_tilt_ij(atoms, c1, c2, fold_atoms) 

426 

427 

428def update_cell_and_positions(atoms, new_cell, op): 

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

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

431 scpos %= 1.0 

432 scpos %= 1.0 

433 

434 atoms.set_cell(new_cell) 

435 atoms.set_scaled_positions(scpos) 

436 

437 

438def niggli_reduce(atoms): 

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

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

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

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

443 

444 References: 

445 

446 Niggli, P. "Krystallographische und strukturtheoretische Grundbegriffe. 

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

448 

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

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

451 

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

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

454 2004, A60, 1-6. 

455 """ 

456 

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

458 new_cell, op = niggli_reduce_cell(atoms.cell) 

459 update_cell_and_positions(atoms, new_cell, op) 

460 

461 

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

463 """Reduce atoms object to canonical lattice. 

464 

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

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

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

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

469 from ase.lattice import identify_lattice 

470 niggli_reduce(atoms) 

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

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

473 

474 

475def sort(atoms, tags=None): 

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

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

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

479 

480 Example: 

481 

482 >>> from ase.build import bulk 

483 >>> # Two unit cells of NaCl: 

484 >>> a = 5.64 

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

486 >>> nacl.get_chemical_symbols() 

487 ['Na', 'Cl', 'Na', 'Cl'] 

488 >>> nacl_sorted = sort(nacl) 

489 >>> nacl_sorted.get_chemical_symbols() 

490 ['Cl', 'Cl', 'Na', 'Na'] 

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

492 True 

493 """ 

494 if tags is None: 

495 tags = atoms.get_chemical_symbols() 

496 else: 

497 tags = list(tags) 

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

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

500 return atoms[indices]