Coverage for ase / build / surface.py: 84.10%
239 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3"""Helper functions for creating the most common surfaces and related tasks.
5The helper functions can create the most common low-index surfaces,
6add vacuum layers and add adsorbates.
8"""
10from math import sqrt
11from operator import itemgetter
13import numpy as np
15from ase.atom import Atom
16from ase.atoms import Atoms
17from ase.data import atomic_numbers, reference_states
18from ase.lattice.cubic import FaceCenteredCubic
21def fcc100(symbol, size, a=None, vacuum=None, orthogonal=True,
22 periodic=False):
23 """FCC(100) surface.
25 Supported special adsorption sites: 'ontop', 'bridge', 'hollow'."""
26 if not orthogonal:
27 raise NotImplementedError("Can't do non-orthogonal cell yet!")
29 return _surface(symbol, 'fcc', '100', size, a, None, vacuum,
30 periodic=periodic,
31 orthogonal=orthogonal)
34def fcc110(symbol, size, a=None, vacuum=None, orthogonal=True,
35 periodic=False):
36 """FCC(110) surface.
38 Supported special adsorption sites: 'ontop', 'longbridge',
39 'shortbridge', 'hollow'."""
40 if not orthogonal:
41 raise NotImplementedError("Can't do non-orthogonal cell yet!")
43 return _surface(symbol, 'fcc', '110', size, a, None, vacuum,
44 periodic=periodic,
45 orthogonal=orthogonal)
48def bcc100(symbol, size, a=None, vacuum=None, orthogonal=True,
49 periodic=False):
50 """BCC(100) surface.
52 Supported special adsorption sites: 'ontop', 'bridge', 'hollow'."""
53 if not orthogonal:
54 raise NotImplementedError("Can't do non-orthogonal cell yet!")
56 return _surface(symbol, 'bcc', '100', size, a, None, vacuum,
57 periodic=periodic,
58 orthogonal=orthogonal)
61def bcc110(symbol, size, a=None, vacuum=None, orthogonal=False,
62 periodic=False):
63 """BCC(110) surface.
65 Supported special adsorption sites: 'ontop', 'longbridge',
66 'shortbridge', 'hollow'.
68 Use *orthogonal=True* to get an orthogonal unit cell - works only
69 for size=(i,j,k) with j even."""
70 return _surface(symbol, 'bcc', '110', size, a, None, vacuum,
71 periodic=periodic,
72 orthogonal=orthogonal)
75def bcc111(symbol, size, a=None, vacuum=None, orthogonal=False,
76 periodic=False):
77 """BCC(111) surface.
79 Supported special adsorption sites: 'ontop'.
81 Use *orthogonal=True* to get an orthogonal unit cell - works only
82 for size=(i,j,k) with j even."""
83 return _surface(symbol, 'bcc', '111', size, a, None, vacuum,
84 periodic=periodic,
85 orthogonal=orthogonal)
88def fcc111(symbol, size, a=None, vacuum=None, orthogonal=False,
89 periodic=False):
90 """FCC(111) surface.
92 Supported special adsorption sites: 'ontop', 'bridge', 'fcc' and 'hcp'.
94 Use *orthogonal=True* to get an orthogonal unit cell - works only
95 for size=(i,j,k) with j even."""
96 return _surface(symbol, 'fcc', '111', size, a, None, vacuum,
97 periodic=periodic,
98 orthogonal=orthogonal)
101def hcp0001(symbol, size, a=None, c=None, vacuum=None, orthogonal=False,
102 periodic=False):
103 """HCP(0001) surface.
105 Supported special adsorption sites: 'ontop', 'bridge', 'fcc' and 'hcp'.
107 Use *orthogonal=True* to get an orthogonal unit cell - works only
108 for size=(i,j,k) with j even."""
109 return _surface(symbol, 'hcp', '0001', size, a, c, vacuum,
110 periodic=periodic,
111 orthogonal=orthogonal)
114def hcp10m10(symbol, size, a=None, c=None, vacuum=None, orthogonal=True,
115 periodic=False):
116 """HCP(10m10) surface.
118 Supported special adsorption sites: 'ontop'.
120 Works only for size=(i,j,k) with j even."""
121 if not orthogonal:
122 raise NotImplementedError("Can't do non-orthogonal cell yet!")
124 return _surface(symbol, 'hcp', '10m10', size, a, c, vacuum,
125 periodic=periodic,
126 orthogonal=orthogonal)
129def diamond100(symbol, size, a=None, vacuum=None, orthogonal=True,
130 periodic=False):
131 """DIAMOND(100) surface.
133 Supported special adsorption sites: 'ontop'."""
134 if not orthogonal:
135 raise NotImplementedError("Can't do non-orthogonal cell yet!")
137 return _surface(symbol, 'diamond', '100', size, a, None, vacuum,
138 periodic=periodic,
139 orthogonal=orthogonal)
142def diamond111(symbol, size, a=None, vacuum=None, orthogonal=False,
143 periodic=False):
144 """DIAMOND(111) surface.
146 Supported special adsorption sites: 'ontop'."""
148 if orthogonal:
149 raise NotImplementedError("Can't do orthogonal cell yet!")
150 return _surface(symbol, 'diamond', '111', size, a, None, vacuum,
151 periodic=periodic,
152 orthogonal=orthogonal)
155def add_adsorbate(slab, adsorbate, height, position=(0, 0), offset=None,
156 mol_index=0):
157 """Add an adsorbate to a surface.
159 This function adds an adsorbate to a slab. If the slab is
160 produced by one of the utility functions in ase.build, it
161 is possible to specify the position of the adsorbate by a keyword
162 (the supported keywords depend on which function was used to
163 create the slab).
165 If the adsorbate is a molecule, the atom indexed by the mol_index
166 optional argument is positioned on top of the adsorption position
167 on the surface, and it is the responsibility of the user to orient
168 the adsorbate in a sensible way.
170 This function can be called multiple times to add more than one
171 adsorbate.
173 Parameters
174 ----------
176 slab: The surface onto which the adsorbate should be added.
178 adsorbate: The adsorbate. Must be one of the following three types:
179 A string containing the chemical symbol for a single atom.
180 An atom object.
181 An atoms object (for a molecular adsorbate).
183 height: Height above the surface.
185 position: The x-y position of the adsorbate, either as a tuple of
186 two numbers or as a keyword (if the surface is produced by one
187 of the functions in ase.build).
189 offset (default: None): Offsets the adsorbate by a number of unit
190 cells. Mostly useful when adding more than one adsorbate.
192 mol_index (default: 0): If the adsorbate is a molecule, index of
193 the atom to be positioned above the location specified by the
194 position argument.
196 Note *position* is given in absolute xy coordinates (or as
197 a keyword), whereas offset is specified in unit cells. This
198 can be used to give the positions in units of the unit cell by
199 using *offset* instead.
201 """
202 info = slab.info.get('adsorbate_info', {})
204 pos = np.array([0.0, 0.0]) # (x, y) part
205 spos = np.array([0.0, 0.0]) # part relative to unit cell
206 if offset is not None:
207 spos += np.asarray(offset, float)
209 if isinstance(position, str):
210 # A site-name:
211 if 'sites' not in info:
212 raise TypeError('If the atoms are not made by an ' +
213 'ase.build function, ' +
214 'position cannot be a name.')
215 if position not in info['sites']:
216 raise TypeError(f'Adsorption site {position} not supported.')
217 spos += info['sites'][position]
218 else:
219 pos += position
221 if 'cell' in info:
222 cell = info['cell']
223 else:
224 cell = slab.get_cell()[:2, :2]
226 pos += np.dot(spos, cell)
228 # Convert the adsorbate to an Atoms object
229 if isinstance(adsorbate, Atoms):
230 ads = adsorbate
231 elif isinstance(adsorbate, Atom):
232 ads = Atoms([adsorbate])
233 else:
234 # Assume it is a string representing a single Atom
235 ads = Atoms([Atom(adsorbate)])
237 # Get the z-coordinate:
238 if 'top layer atom index' in info:
239 a = info['top layer atom index']
240 else:
241 a = slab.positions[:, 2].argmax()
242 if 'adsorbate_info' not in slab.info:
243 slab.info['adsorbate_info'] = {}
244 slab.info['adsorbate_info']['top layer atom index'] = a
245 z = slab.positions[a, 2] + height
247 # Move adsorbate into position
248 ads.translate([pos[0], pos[1], z] - ads.positions[mol_index])
250 # Attach the adsorbate
251 slab.extend(ads)
254def add_vacuum(atoms, vacuum):
255 """Add vacuum layer to the atoms.
257 Parameters
258 ----------
260 atoms: Atoms object
261 Most likely created by one of the surface functions.
262 vacuum: float
263 The thickness of the vacuum layer (in Angstrom).
264 """
265 uc = atoms.get_cell()
266 normal = np.cross(uc[0], uc[1])
267 costheta = np.dot(normal, uc[2]) / np.sqrt(np.dot(normal, normal) *
268 np.dot(uc[2], uc[2]))
269 length = np.sqrt(np.dot(uc[2], uc[2]))
270 newlength = length + vacuum / costheta
271 uc[2] *= newlength / length
272 atoms.set_cell(uc)
275def create_tags(size) -> np.ndarray:
276 """ Function to create layer tags. """
277 # tag atoms by layer
278 # create blocks of descending integers of length size[0]*size[1]
279 return np.arange(size[2], 0, -1).repeat(size[0] * size[1])
282def _surface(symbol, structure, face, size, a, c, vacuum, periodic,
283 orthogonal=True):
284 """Function to build often used surfaces.
286 Don't call this function directly - use fcc100, fcc110, bcc111, ..."""
288 Z = atomic_numbers[symbol]
290 if a is None:
291 sym = reference_states[Z]['symmetry']
292 if sym != structure:
293 raise ValueError(
294 f"Can't guess lattice constant for {structure}-{symbol}!")
295 a = reference_states[Z]['a']
297 if structure == 'hcp' and c is None:
298 if reference_states[Z]['symmetry'] == 'hcp':
299 c = reference_states[Z]['c/a'] * a
300 else:
301 c = sqrt(8 / 3.0) * a
303 positions = np.empty((size[2], size[1], size[0], 3))
304 positions[..., 0] = np.arange(size[0]).reshape((1, 1, -1))
305 positions[..., 1] = np.arange(size[1]).reshape((1, -1, 1))
306 positions[..., 2] = np.arange(size[2]).reshape((-1, 1, 1))
308 numbers = np.ones(size[0] * size[1] * size[2], int) * Z
310 slab = Atoms(numbers,
311 tags=create_tags(size),
312 pbc=(True, True, periodic),
313 cell=size)
315 surface_cell = None
316 sites = {'ontop': (0, 0)}
317 surf = structure + face
318 if surf == 'fcc100':
319 cell = (sqrt(0.5), sqrt(0.5), 0.5)
320 positions[-2::-2, ..., :2] += 0.5
321 sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)})
322 elif surf == 'diamond100':
323 cell = (sqrt(0.5), sqrt(0.5), 0.5 / 2)
324 positions[-4::-4, ..., :2] += (0.5, 0.5)
325 positions[-3::-4, ..., :2] += (0.0, 0.5)
326 positions[-2::-4, ..., :2] += (0.0, 0.0)
327 positions[-1::-4, ..., :2] += (0.5, 0.0)
328 elif surf == 'fcc110':
329 cell = (1.0, sqrt(0.5), sqrt(0.125))
330 positions[-2::-2, ..., :2] += 0.5
331 sites.update({'hollow': (0.5, 0.5), 'longbridge': (0.5, 0),
332 'shortbridge': (0, 0.5)})
333 elif surf == 'bcc100':
334 cell = (1.0, 1.0, 0.5)
335 positions[-2::-2, ..., :2] += 0.5
336 sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)})
337 else:
338 if orthogonal and size[1] % 2 == 1:
339 raise ValueError(("Can't make orthorhombic cell with size=%r. " %
340 (tuple(size),)) +
341 'Second number in size must be even.')
342 if surf == 'fcc111':
343 cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3))
344 if orthogonal:
345 positions[-1::-3, 1::2, :, 0] += 0.5
346 positions[-2::-3, 1::2, :, 0] += 0.5
347 positions[-3::-3, 1::2, :, 0] -= 0.5
348 positions[-2::-3, ..., :2] += (0.0, 2.0 / 3)
349 positions[-3::-3, ..., :2] += (0.5, 1.0 / 3)
350 else:
351 positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3)
352 positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3)
353 sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3),
354 'hcp': (2.0 / 3, 2.0 / 3)})
355 elif surf == 'diamond111':
356 cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3) / 2)
357 assert not orthogonal
358 positions[-1::-6, ..., :3] += (0.0, 0.0, 0.5)
359 positions[-2::-6, ..., :2] += (0.0, 0.0)
360 positions[-3::-6, ..., :3] += (-1.0 / 3, 2.0 / 3, 0.5)
361 positions[-4::-6, ..., :2] += (-1.0 / 3, 2.0 / 3)
362 positions[-5::-6, ..., :3] += (1.0 / 3, 1.0 / 3, 0.5)
363 positions[-6::-6, ..., :2] += (1.0 / 3, 1.0 / 3)
364 elif surf == 'hcp0001':
365 cell = (1.0, sqrt(0.75), 0.5 * c / a)
366 if orthogonal:
367 positions[:, 1::2, :, 0] += 0.5
368 positions[-2::-2, ..., :2] += (0.0, 2.0 / 3)
369 else:
370 positions[-2::-2, ..., :2] += (-1.0 / 3, 2.0 / 3)
371 sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3),
372 'hcp': (2.0 / 3, 2.0 / 3)})
373 elif surf == 'hcp10m10':
374 cell = (1.0, 0.5 * c / a, sqrt(0.75))
375 assert orthogonal
376 positions[-2::-2, ..., 0] += 0.5
377 positions[:, ::2, :, 2] += 2.0 / 3
378 elif surf == 'bcc110':
379 cell = (1.0, sqrt(0.5), sqrt(0.5))
380 if orthogonal:
381 positions[:, 1::2, :, 0] += 0.5
382 positions[-2::-2, ..., :2] += (0.0, 1.0)
383 else:
384 positions[-2::-2, ..., :2] += (-0.5, 1.0)
385 sites.update({'shortbridge': (0, 0.5),
386 'longbridge': (0.5, 0),
387 'hollow': (0.375, 0.25)})
388 elif surf == 'bcc111':
389 cell = (sqrt(2), sqrt(1.5), sqrt(3) / 6)
390 if orthogonal:
391 positions[-1::-3, 1::2, :, 0] += 0.5
392 positions[-2::-3, 1::2, :, 0] += 0.5
393 positions[-3::-3, 1::2, :, 0] -= 0.5
394 positions[-2::-3, ..., :2] += (0.0, 2.0 / 3)
395 positions[-3::-3, ..., :2] += (0.5, 1.0 / 3)
396 else:
397 positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3)
398 positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3)
399 sites.update({'hollow': (1.0 / 3, 1.0 / 3)})
400 else:
401 2 / 0
403 surface_cell = a * np.array([(cell[0], 0),
404 (cell[0] / 2, cell[1])])
405 if not orthogonal:
406 cell = np.array([(cell[0], 0, 0),
407 (cell[0] / 2, cell[1], 0),
408 (0, 0, cell[2])])
410 if surface_cell is None:
411 surface_cell = a * np.diag(cell[:2])
413 if isinstance(cell, tuple):
414 cell = np.diag(cell)
416 slab.set_positions(positions.reshape((-1, 3)))
417 slab.set_cell([a * v * n for v, n in zip(cell, size)], scale_atoms=True)
419 if not periodic:
420 slab.cell[2] = 0.0
422 if vacuum is not None:
423 slab.center(vacuum, axis=2)
425 if 'adsorbate_info' not in slab.info:
426 slab.info.update({'adsorbate_info': {}})
428 slab.info['adsorbate_info']['cell'] = surface_cell
429 slab.info['adsorbate_info']['sites'] = sites
430 return slab
433def fcc211(symbol, size, a=None, vacuum=None, orthogonal=True):
434 """FCC(211) surface.
436 Does not currently support special adsorption sites.
438 Currently only implemented for *orthogonal=True* with size specified
439 as (i, j, k), where i, j, and k are number of atoms in each direction.
440 i must be divisible by 3 to accommodate the step width.
441 """
442 if not orthogonal:
443 raise NotImplementedError('Only implemented for orthogonal '
444 'unit cells.')
445 if size[0] % 3 != 0:
446 raise NotImplementedError('First dimension of size must be '
447 'divisible by 3.')
448 atoms = FaceCenteredCubic(symbol,
449 directions=[[1, -1, -1],
450 [0, 2, -2],
451 [2, 1, 1]],
452 miller=(None, None, (2, 1, 1)),
453 latticeconstant=a,
454 size=(1, 1, 1),
455 pbc=True)
456 z = (size[2] + 1) // 2
457 atoms = atoms.repeat((size[0] // 3, size[1], z))
458 if size[2] % 2: # Odd: remove bottom layer and shrink cell.
459 remove_list = [atom.index for atom in atoms
460 if atom.z < atoms[1].z]
461 del atoms[remove_list]
462 dz = atoms[0].z
463 atoms.translate((0., 0., -dz))
464 atoms.cell[2][2] -= dz
466 atoms.cell[2] = 0.0
467 atoms.pbc[2] = False
468 if vacuum:
469 atoms.center(vacuum, axis=2)
471 # Renumber systematically from top down.
472 orders = [(atom.index, round(atom.x, 3), round(atom.y, 3),
473 -round(atom.z, 3), atom.index) for atom in atoms]
474 orders.sort(key=itemgetter(3, 1, 2))
475 newatoms = atoms.copy()
476 for index, order in enumerate(orders):
477 newatoms[index].position = atoms[order[0]].position.copy()
479 # Add empty 'sites' dictionary for consistency with other functions
480 newatoms.info['adsorbate_info'] = {'sites': {}}
481 return newatoms
484def mx2(formula='MoS2', kind='2H', a=3.18, thickness=3.19,
485 size=(1, 1, 1), vacuum=None):
486 """Create three-layer 2D materials with hexagonal structure.
488 This can be used for e.g. metal dichalcogenides :mol:`MX_2` 2D structures
489 such as :mol:`MoS_2`.
491 https://en.wikipedia.org/wiki/Transition_metal_dichalcogenide_monolayers
493 Parameters
494 ----------
495 kind : {'2H', '1T'}, default: '2H'
497 - '2H': mirror-plane symmetry
498 - '1T': inversion symmetry
499 """
500 if kind == '2H':
501 basis = [(0, 0, 0),
502 (2 / 3, 1 / 3, 0.5 * thickness),
503 (2 / 3, 1 / 3, -0.5 * thickness)]
504 elif kind == '1T':
505 basis = [(0, 0, 0),
506 (2 / 3, 1 / 3, 0.5 * thickness),
507 (1 / 3, 2 / 3, -0.5 * thickness)]
508 else:
509 raise ValueError('Structure not recognized:', kind)
511 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
513 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
514 atoms.set_scaled_positions(basis)
515 if vacuum is not None:
516 atoms.center(vacuum, axis=2)
517 atoms = atoms.repeat(size)
518 return atoms
521def graphene(formula='C2', a=2.460, thickness=0.0,
522 size=(1, 1, 1), vacuum=None):
523 """Create a graphene monolayer structure.
525 Parameters
526 ----------
527 thickness : float, default: 0.0
528 Thickness of the layer; maybe for a buckled structure like silicene.
529 """
530 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
531 basis = [[0, 0, -0.5 * thickness], [2 / 3, 1 / 3, 0.5 * thickness]]
532 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
533 atoms.set_scaled_positions(basis)
534 if vacuum is not None:
535 atoms.center(vacuum, axis=2)
536 atoms = atoms.repeat(size)
537 return atoms
540def _all_surface_functions():
541 # Convenient for debugging.
542 d = {
543 func.__name__: func
544 for func in [
545 fcc100,
546 fcc110,
547 bcc100,
548 bcc110,
549 bcc111,
550 fcc111,
551 hcp0001,
552 hcp10m10,
553 diamond100,
554 diamond111,
555 fcc111,
556 mx2,
557 graphene,
558 ]
559 }
560 return d