 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np

2from ase.build.niggli import niggli_reduce_cell

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*.

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

19 Parameters:

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*.

74 Example:

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

107 if isinstance(origo, int):

108 origo = atoms.get_scaled_positions()[origo]

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

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)

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

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)

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

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

155 at.cell *= levels[nlayers]

156 return at[tags < nlayers]

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)

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)

176 # Mask out atoms outside new cell

177 stol = 0.1 * tolerance # scaled tolerance, XXX

178 maskcell = atoms.cell * extend

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

182 return atoms

185class IncompatibleCellError(ValueError):

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

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

188 pass

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.

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*.

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.

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.

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'.

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

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

225 structure.

227 Example:

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

258 for atoms in [atoms1, atoms2]:

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

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

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

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

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

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

295 if distance is not None:

296 from scipy.optimize import fmin

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

305 def func(x):

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

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

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, x

318 atoms1.translate(t1)

319 atoms2.translate(t2)

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

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

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

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

325 atoms1.extend(atoms2)

327 if reorder:

328 atoms1 = sort(atoms1)

330 if output_strained:

331 return atoms1, atoms1_strained, atoms2_strained

332 else:

333 return atoms1

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*.

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

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

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

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

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.

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

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

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

383 if rotate_cell:

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

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

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

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

391 that is kept fixed.

392 """

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

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]

406 # sanity check

407 def volume(cell):

408 return np.abs(np.dot(cell, np.cross(cell, cell)))

409 V = volume(cell_cc)

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

412 atoms.set_cell(cell_cc)

414 if fold_atoms:

415 atoms.wrap()

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

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)

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

434 atoms.set_cell(new_cell)

435 atoms.set_scaled_positions(scpos)

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.

444 References:

446 Niggli, P. "Krystallographische und strukturtheoretische Grundbegriffe.

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

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

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

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

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)

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

463 """Reduce atoms object to canonical lattice.

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

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.

480 Example:

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]