1"""Helper functions for creating the most common surfaces and related tasks. 

2 

3The helper functions can create the most common low-index surfaces, 

4add vacuum layers and add adsorbates. 

5 

6""" 

7 

8from math import sqrt 

9from operator import itemgetter 

10 

11import numpy as np 

12 

13from ase.atom import Atom 

14from ase.atoms import Atoms 

15from ase.data import reference_states, atomic_numbers 

16from ase.lattice.cubic import FaceCenteredCubic 

17 

18 

19def fcc100(symbol, size, a=None, vacuum=None, orthogonal=True, 

20 periodic=False): 

21 """FCC(100) surface. 

22 

23 Supported special adsorption sites: 'ontop', 'bridge', 'hollow'.""" 

24 if not orthogonal: 

25 raise NotImplementedError("Can't do non-orthogonal cell yet!") 

26 

27 return _surface(symbol, 'fcc', '100', size, a, None, vacuum, 

28 periodic=periodic, 

29 orthogonal=orthogonal) 

30 

31 

32def fcc110(symbol, size, a=None, vacuum=None, orthogonal=True, 

33 periodic=False): 

34 """FCC(110) surface. 

35 

36 Supported special adsorption sites: 'ontop', 'longbridge', 

37 'shortbridge', 'hollow'.""" 

38 if not orthogonal: 

39 raise NotImplementedError("Can't do non-orthogonal cell yet!") 

40 

41 return _surface(symbol, 'fcc', '110', size, a, None, vacuum, 

42 periodic=periodic, 

43 orthogonal=orthogonal) 

44 

45 

46def bcc100(symbol, size, a=None, vacuum=None, orthogonal=True, 

47 periodic=False): 

48 """BCC(100) surface. 

49 

50 Supported special adsorption sites: 'ontop', 'bridge', 'hollow'.""" 

51 if not orthogonal: 

52 raise NotImplementedError("Can't do non-orthogonal cell yet!") 

53 

54 return _surface(symbol, 'bcc', '100', size, a, None, vacuum, 

55 periodic=periodic, 

56 orthogonal=orthogonal) 

57 

58 

59def bcc110(symbol, size, a=None, vacuum=None, orthogonal=False, 

60 periodic=False): 

61 """BCC(110) surface. 

62 

63 Supported special adsorption sites: 'ontop', 'longbridge', 

64 'shortbridge', 'hollow'. 

65 

66 Use *orthogonal=True* to get an orthogonal unit cell - works only 

67 for size=(i,j,k) with j even.""" 

68 return _surface(symbol, 'bcc', '110', size, a, None, vacuum, 

69 periodic=periodic, 

70 orthogonal=orthogonal) 

71 

72 

73def bcc111(symbol, size, a=None, vacuum=None, orthogonal=False, 

74 periodic=False): 

75 """BCC(111) surface. 

76 

77 Supported special adsorption sites: 'ontop'. 

78 

79 Use *orthogonal=True* to get an orthogonal unit cell - works only 

80 for size=(i,j,k) with j even.""" 

81 return _surface(symbol, 'bcc', '111', size, a, None, vacuum, 

82 periodic=periodic, 

83 orthogonal=orthogonal) 

84 

85 

86def fcc111(symbol, size, a=None, vacuum=None, orthogonal=False, 

87 periodic=False): 

88 """FCC(111) surface. 

89 

90 Supported special adsorption sites: 'ontop', 'bridge', 'fcc' and 'hcp'. 

91 

92 Use *orthogonal=True* to get an orthogonal unit cell - works only 

93 for size=(i,j,k) with j even.""" 

94 return _surface(symbol, 'fcc', '111', size, a, None, vacuum, 

95 periodic=periodic, 

96 orthogonal=orthogonal) 

97 

98 

99def hcp0001(symbol, size, a=None, c=None, vacuum=None, orthogonal=False, 

100 periodic=False): 

101 """HCP(0001) surface. 

102 

103 Supported special adsorption sites: 'ontop', 'bridge', 'fcc' and 'hcp'. 

104 

105 Use *orthogonal=True* to get an orthogonal unit cell - works only 

106 for size=(i,j,k) with j even.""" 

107 return _surface(symbol, 'hcp', '0001', size, a, c, vacuum, 

108 periodic=periodic, 

109 orthogonal=orthogonal) 

110 

111 

112def hcp10m10(symbol, size, a=None, c=None, vacuum=None, orthogonal=True, 

113 periodic=False): 

114 """HCP(10m10) surface. 

115 

116 Supported special adsorption sites: 'ontop'. 

117 

118 Works only for size=(i,j,k) with j even.""" 

119 if not orthogonal: 

120 raise NotImplementedError("Can't do non-orthogonal cell yet!") 

121 

122 return _surface(symbol, 'hcp', '10m10', size, a, c, vacuum, 

123 periodic=periodic, 

124 orthogonal=orthogonal) 

125 

126 

127def diamond100(symbol, size, a=None, vacuum=None, orthogonal=True, 

128 periodic=False): 

129 """DIAMOND(100) surface. 

130 

131 Supported special adsorption sites: 'ontop'.""" 

132 if not orthogonal: 

133 raise NotImplementedError("Can't do non-orthogonal cell yet!") 

134 

135 return _surface(symbol, 'diamond', '100', size, a, None, vacuum, 

136 periodic=periodic, 

137 orthogonal=orthogonal) 

138 

139 

140def diamond111(symbol, size, a=None, vacuum=None, orthogonal=False, 

141 periodic=False): 

142 """DIAMOND(111) surface. 

143 

144 Supported special adsorption sites: 'ontop'.""" 

145 

146 if orthogonal: 

147 raise NotImplementedError("Can't do orthogonal cell yet!") 

148 return _surface(symbol, 'diamond', '111', size, a, None, vacuum, 

149 periodic=periodic, 

150 orthogonal=orthogonal) 

151 

152 

153def add_adsorbate(slab, adsorbate, height, position=(0, 0), offset=None, 

154 mol_index=0): 

155 """Add an adsorbate to a surface. 

156 

157 This function adds an adsorbate to a slab. If the slab is 

158 produced by one of the utility functions in ase.build, it 

159 is possible to specify the position of the adsorbate by a keyword 

160 (the supported keywords depend on which function was used to 

161 create the slab). 

162 

163 If the adsorbate is a molecule, the atom indexed by the mol_index 

164 optional argument is positioned on top of the adsorption position 

165 on the surface, and it is the responsibility of the user to orient 

166 the adsorbate in a sensible way. 

167 

168 This function can be called multiple times to add more than one 

169 adsorbate. 

170 

171 Parameters: 

172 

173 slab: The surface onto which the adsorbate should be added. 

174 

175 adsorbate: The adsorbate. Must be one of the following three types: 

176 A string containing the chemical symbol for a single atom. 

177 An atom object. 

178 An atoms object (for a molecular adsorbate). 

179 

180 height: Height above the surface. 

181 

182 position: The x-y position of the adsorbate, either as a tuple of 

183 two numbers or as a keyword (if the surface is produced by one 

184 of the functions in ase.build). 

185 

186 offset (default: None): Offsets the adsorbate by a number of unit 

187 cells. Mostly useful when adding more than one adsorbate. 

188 

189 mol_index (default: 0): If the adsorbate is a molecule, index of 

190 the atom to be positioned above the location specified by the 

191 position argument. 

192 

193 Note *position* is given in absolute xy coordinates (or as 

194 a keyword), whereas offset is specified in unit cells. This 

195 can be used to give the positions in units of the unit cell by 

196 using *offset* instead. 

197 

198 """ 

199 info = slab.info.get('adsorbate_info', {}) 

200 

201 pos = np.array([0.0, 0.0]) # (x, y) part 

202 spos = np.array([0.0, 0.0]) # part relative to unit cell 

203 if offset is not None: 

204 spos += np.asarray(offset, float) 

205 

206 if isinstance(position, str): 

207 # A site-name: 

208 if 'sites' not in info: 

209 raise TypeError('If the atoms are not made by an ' + 

210 'ase.build function, ' + 

211 'position cannot be a name.') 

212 if position not in info['sites']: 

213 raise TypeError('Adsorption site %s not supported.' % position) 

214 spos += info['sites'][position] 

215 else: 

216 pos += position 

217 

218 if 'cell' in info: 

219 cell = info['cell'] 

220 else: 

221 cell = slab.get_cell()[:2, :2] 

222 

223 pos += np.dot(spos, cell) 

224 

225 # Convert the adsorbate to an Atoms object 

226 if isinstance(adsorbate, Atoms): 

227 ads = adsorbate 

228 elif isinstance(adsorbate, Atom): 

229 ads = Atoms([adsorbate]) 

230 else: 

231 # Assume it is a string representing a single Atom 

232 ads = Atoms([Atom(adsorbate)]) 

233 

234 # Get the z-coordinate: 

235 if 'top layer atom index' in info: 

236 a = info['top layer atom index'] 

237 else: 

238 a = slab.positions[:, 2].argmax() 

239 if 'adsorbate_info' not in slab.info: 

240 slab.info['adsorbate_info'] = {} 

241 slab.info['adsorbate_info']['top layer atom index'] = a 

242 z = slab.positions[a, 2] + height 

243 

244 # Move adsorbate into position 

245 ads.translate([pos[0], pos[1], z] - ads.positions[mol_index]) 

246 

247 # Attach the adsorbate 

248 slab.extend(ads) 

249 

250 

251def add_vacuum(atoms, vacuum): 

252 """Add vacuum layer to the atoms. 

253 

254 Parameters: 

255 

256 atoms: Atoms object 

257 Most likely created by one of the surface functions. 

258 vacuum: float 

259 The thickness of the vacuum layer (in Angstrom). 

260 """ 

261 uc = atoms.get_cell() 

262 normal = np.cross(uc[0], uc[1]) 

263 costheta = np.dot(normal, uc[2]) / np.sqrt(np.dot(normal, normal) * 

264 np.dot(uc[2], uc[2])) 

265 length = np.sqrt(np.dot(uc[2], uc[2])) 

266 newlength = length + vacuum / costheta 

267 uc[2] *= newlength / length 

268 atoms.set_cell(uc) 

269 

270 

271def _surface(symbol, structure, face, size, a, c, vacuum, periodic, 

272 orthogonal=True): 

273 """Function to build often used surfaces. 

274 

275 Don't call this function directly - use fcc100, fcc110, bcc111, ...""" 

276 

277 Z = atomic_numbers[symbol] 

278 

279 if a is None: 

280 sym = reference_states[Z]['symmetry'] 

281 if sym != structure: 

282 raise ValueError("Can't guess lattice constant for %s-%s!" % 

283 (structure, symbol)) 

284 a = reference_states[Z]['a'] 

285 

286 if structure == 'hcp' and c is None: 

287 if reference_states[Z]['symmetry'] == 'hcp': 

288 c = reference_states[Z]['c/a'] * a 

289 else: 

290 c = sqrt(8 / 3.0) * a 

291 

292 positions = np.empty((size[2], size[1], size[0], 3)) 

293 positions[..., 0] = np.arange(size[0]).reshape((1, 1, -1)) 

294 positions[..., 1] = np.arange(size[1]).reshape((1, -1, 1)) 

295 positions[..., 2] = np.arange(size[2]).reshape((-1, 1, 1)) 

296 

297 numbers = np.ones(size[0] * size[1] * size[2], int) * Z 

298 

299 tags = np.empty((size[2], size[1], size[0]), int) 

300 tags[:] = np.arange(size[2], 0, -1).reshape((-1, 1, 1)) 

301 

302 slab = Atoms(numbers, 

303 tags=tags.ravel(), 

304 pbc=(True, True, periodic), 

305 cell=size) 

306 

307 surface_cell = None 

308 sites = {'ontop': (0, 0)} 

309 surf = structure + face 

310 if surf == 'fcc100': 

311 cell = (sqrt(0.5), sqrt(0.5), 0.5) 

312 positions[-2::-2, ..., :2] += 0.5 

313 sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)}) 

314 elif surf == 'diamond100': 

315 cell = (sqrt(0.5), sqrt(0.5), 0.5 / 2) 

316 positions[-4::-4, ..., :2] += (0.5, 0.5) 

317 positions[-3::-4, ..., :2] += (0.0, 0.5) 

318 positions[-2::-4, ..., :2] += (0.0, 0.0) 

319 positions[-1::-4, ..., :2] += (0.5, 0.0) 

320 elif surf == 'fcc110': 

321 cell = (1.0, sqrt(0.5), sqrt(0.125)) 

322 positions[-2::-2, ..., :2] += 0.5 

323 sites.update({'hollow': (0.5, 0.5), 'longbridge': (0.5, 0), 

324 'shortbridge': (0, 0.5)}) 

325 elif surf == 'bcc100': 

326 cell = (1.0, 1.0, 0.5) 

327 positions[-2::-2, ..., :2] += 0.5 

328 sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)}) 

329 else: 

330 if orthogonal and size[1] % 2 == 1: 

331 raise ValueError(("Can't make orthorhombic cell with size=%r. " % 

332 (tuple(size),)) + 

333 'Second number in size must be even.') 

334 if surf == 'fcc111': 

335 cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3)) 

336 if orthogonal: 

337 positions[-1::-3, 1::2, :, 0] += 0.5 

338 positions[-2::-3, 1::2, :, 0] += 0.5 

339 positions[-3::-3, 1::2, :, 0] -= 0.5 

340 positions[-2::-3, ..., :2] += (0.0, 2.0 / 3) 

341 positions[-3::-3, ..., :2] += (0.5, 1.0 / 3) 

342 else: 

343 positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3) 

344 positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3) 

345 sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3), 

346 'hcp': (2.0 / 3, 2.0 / 3)}) 

347 elif surf == 'diamond111': 

348 cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3) / 2) 

349 assert not orthogonal 

350 positions[-1::-6, ..., :3] += (0.0, 0.0, 0.5) 

351 positions[-2::-6, ..., :2] += (0.0, 0.0) 

352 positions[-3::-6, ..., :3] += (-1.0 / 3, 2.0 / 3, 0.5) 

353 positions[-4::-6, ..., :2] += (-1.0 / 3, 2.0 / 3) 

354 positions[-5::-6, ..., :3] += (1.0 / 3, 1.0 / 3, 0.5) 

355 positions[-6::-6, ..., :2] += (1.0 / 3, 1.0 / 3) 

356 elif surf == 'hcp0001': 

357 cell = (1.0, sqrt(0.75), 0.5 * c / a) 

358 if orthogonal: 

359 positions[:, 1::2, :, 0] += 0.5 

360 positions[-2::-2, ..., :2] += (0.0, 2.0 / 3) 

361 else: 

362 positions[-2::-2, ..., :2] += (-1.0 / 3, 2.0 / 3) 

363 sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3), 

364 'hcp': (2.0 / 3, 2.0 / 3)}) 

365 elif surf == 'hcp10m10': 

366 cell = (1.0, 0.5 * c / a, sqrt(0.75)) 

367 assert orthogonal 

368 positions[-2::-2, ..., 0] += 0.5 

369 positions[:, ::2, :, 2] += 2.0 / 3 

370 elif surf == 'bcc110': 

371 cell = (1.0, sqrt(0.5), sqrt(0.5)) 

372 if orthogonal: 

373 positions[:, 1::2, :, 0] += 0.5 

374 positions[-2::-2, ..., :2] += (0.0, 1.0) 

375 else: 

376 positions[-2::-2, ..., :2] += (-0.5, 1.0) 

377 sites.update({'shortbridge': (0, 0.5), 

378 'longbridge': (0.5, 0), 

379 'hollow': (0.375, 0.25)}) 

380 elif surf == 'bcc111': 

381 cell = (sqrt(2), sqrt(1.5), sqrt(3) / 6) 

382 if orthogonal: 

383 positions[-1::-3, 1::2, :, 0] += 0.5 

384 positions[-2::-3, 1::2, :, 0] += 0.5 

385 positions[-3::-3, 1::2, :, 0] -= 0.5 

386 positions[-2::-3, ..., :2] += (0.0, 2.0 / 3) 

387 positions[-3::-3, ..., :2] += (0.5, 1.0 / 3) 

388 else: 

389 positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3) 

390 positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3) 

391 sites.update({'hollow': (1.0 / 3, 1.0 / 3)}) 

392 else: 

393 2 / 0 

394 

395 surface_cell = a * np.array([(cell[0], 0), 

396 (cell[0] / 2, cell[1])]) 

397 if not orthogonal: 

398 cell = np.array([(cell[0], 0, 0), 

399 (cell[0] / 2, cell[1], 0), 

400 (0, 0, cell[2])]) 

401 

402 if surface_cell is None: 

403 surface_cell = a * np.diag(cell[:2]) 

404 

405 if isinstance(cell, tuple): 

406 cell = np.diag(cell) 

407 

408 slab.set_positions(positions.reshape((-1, 3))) 

409 slab.set_cell([a * v * n for v, n in zip(cell, size)], scale_atoms=True) 

410 

411 if not periodic: 

412 slab.cell[2] = 0.0 

413 

414 if vacuum is not None: 

415 slab.center(vacuum, axis=2) 

416 

417 if 'adsorbate_info' not in slab.info: 

418 slab.info.update({'adsorbate_info': {}}) 

419 

420 slab.info['adsorbate_info']['cell'] = surface_cell 

421 slab.info['adsorbate_info']['sites'] = sites 

422 return slab 

423 

424 

425def fcc211(symbol, size, a=None, vacuum=None, orthogonal=True): 

426 """FCC(211) surface. 

427 

428 Does not currently support special adsorption sites. 

429 

430 Currently only implemented for *orthogonal=True* with size specified 

431 as (i, j, k), where i, j, and k are number of atoms in each direction. 

432 i must be divisible by 3 to accommodate the step width. 

433 """ 

434 if not orthogonal: 

435 raise NotImplementedError('Only implemented for orthogonal ' 

436 'unit cells.') 

437 if size[0] % 3 != 0: 

438 raise NotImplementedError('First dimension of size must be ' 

439 'divisible by 3.') 

440 atoms = FaceCenteredCubic(symbol, 

441 directions=[[1, -1, -1], 

442 [0, 2, -2], 

443 [2, 1, 1]], 

444 miller=(None, None, (2, 1, 1)), 

445 latticeconstant=a, 

446 size=(1, 1, 1), 

447 pbc=True) 

448 z = (size[2] + 1) // 2 

449 atoms = atoms.repeat((size[0] // 3, size[1], z)) 

450 if size[2] % 2: # Odd: remove bottom layer and shrink cell. 

451 remove_list = [atom.index for atom in atoms 

452 if atom.z < atoms[1].z] 

453 del atoms[remove_list] 

454 dz = atoms[0].z 

455 atoms.translate((0., 0., -dz)) 

456 atoms.cell[2][2] -= dz 

457 

458 atoms.cell[2] = 0.0 

459 atoms.pbc[2] = False 

460 if vacuum: 

461 atoms.center(vacuum, axis=2) 

462 

463 # Renumber systematically from top down. 

464 orders = [(atom.index, round(atom.x, 3), round(atom.y, 3), 

465 -round(atom.z, 3), atom.index) for atom in atoms] 

466 orders.sort(key=itemgetter(3, 1, 2)) 

467 newatoms = atoms.copy() 

468 for index, order in enumerate(orders): 

469 newatoms[index].position = atoms[order[0]].position.copy() 

470 

471 # Add empty 'sites' dictionary for consistency with other functions 

472 newatoms.info['adsorbate_info'] = {'sites': {}} 

473 return newatoms 

474 

475 

476def mx2(formula='MoS2', kind='2H', a=3.18, thickness=3.19, 

477 size=(1, 1, 1), vacuum=None): 

478 """Create three-layer 2D materials with hexagonal structure. 

479 

480 For metal dichalcogenites, etc. 

481 

482 The kind argument accepts '2H', which gives a mirror plane symmetry 

483 and '1T', which gives an inversion symmetry.""" 

484 

485 if kind == '2H': 

486 basis = [(0, 0, 0), 

487 (2 / 3, 1 / 3, 0.5 * thickness), 

488 (2 / 3, 1 / 3, -0.5 * thickness)] 

489 elif kind == '1T': 

490 basis = [(0, 0, 0), 

491 (2 / 3, 1 / 3, 0.5 * thickness), 

492 (1 / 3, 2 / 3, -0.5 * thickness)] 

493 else: 

494 raise ValueError('Structure not recognized:', kind) 

495 

496 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]] 

497 

498 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0)) 

499 atoms.set_scaled_positions(basis) 

500 if vacuum is not None: 

501 atoms.center(vacuum, axis=2) 

502 atoms = atoms.repeat(size) 

503 return atoms 

504 

505 

506def graphene(formula='C2', a=2.460, size=(1, 1, 1), vacuum=None): 

507 """Create a graphene monolayer structure.""" 

508 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]] 

509 basis = [[0, 0, 0], [2 / 3, 1 / 3, 0]] 

510 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0)) 

511 atoms.set_scaled_positions(basis) 

512 if vacuum is not None: 

513 atoms.center(vacuum, axis=2) 

514 atoms = atoms.repeat(size) 

515 return atoms 

516 

517 

518def _all_surface_functions(): 

519 # Convenient for debugging. 

520 d = {} 

521 for func in [fcc100, fcc110, bcc100, bcc110, bcc111, fcc111, hcp0001, 

522 hcp10m10, diamond100, diamond111, fcc111, mx2, graphene]: 

523 d[func.__name__] = func 

524 return d