Coverage for /builds/ase/ase/ase/build/surface.py: 84.10%
239 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +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:
175 slab: The surface onto which the adsorbate should be added.
177 adsorbate: The adsorbate. Must be one of the following three types:
178 A string containing the chemical symbol for a single atom.
179 An atom object.
180 An atoms object (for a molecular adsorbate).
182 height: Height above the surface.
184 position: The x-y position of the adsorbate, either as a tuple of
185 two numbers or as a keyword (if the surface is produced by one
186 of the functions in ase.build).
188 offset (default: None): Offsets the adsorbate by a number of unit
189 cells. Mostly useful when adding more than one adsorbate.
191 mol_index (default: 0): If the adsorbate is a molecule, index of
192 the atom to be positioned above the location specified by the
193 position argument.
195 Note *position* is given in absolute xy coordinates (or as
196 a keyword), whereas offset is specified in unit cells. This
197 can be used to give the positions in units of the unit cell by
198 using *offset* instead.
200 """
201 info = slab.info.get('adsorbate_info', {})
203 pos = np.array([0.0, 0.0]) # (x, y) part
204 spos = np.array([0.0, 0.0]) # part relative to unit cell
205 if offset is not None:
206 spos += np.asarray(offset, float)
208 if isinstance(position, str):
209 # A site-name:
210 if 'sites' not in info:
211 raise TypeError('If the atoms are not made by an ' +
212 'ase.build function, ' +
213 'position cannot be a name.')
214 if position not in info['sites']:
215 raise TypeError(f'Adsorption site {position} not supported.')
216 spos += info['sites'][position]
217 else:
218 pos += position
220 if 'cell' in info:
221 cell = info['cell']
222 else:
223 cell = slab.get_cell()[:2, :2]
225 pos += np.dot(spos, cell)
227 # Convert the adsorbate to an Atoms object
228 if isinstance(adsorbate, Atoms):
229 ads = adsorbate
230 elif isinstance(adsorbate, Atom):
231 ads = Atoms([adsorbate])
232 else:
233 # Assume it is a string representing a single Atom
234 ads = Atoms([Atom(adsorbate)])
236 # Get the z-coordinate:
237 if 'top layer atom index' in info:
238 a = info['top layer atom index']
239 else:
240 a = slab.positions[:, 2].argmax()
241 if 'adsorbate_info' not in slab.info:
242 slab.info['adsorbate_info'] = {}
243 slab.info['adsorbate_info']['top layer atom index'] = a
244 z = slab.positions[a, 2] + height
246 # Move adsorbate into position
247 ads.translate([pos[0], pos[1], z] - ads.positions[mol_index])
249 # Attach the adsorbate
250 slab.extend(ads)
253def add_vacuum(atoms, vacuum):
254 """Add vacuum layer to the atoms.
256 Parameters:
258 atoms: Atoms object
259 Most likely created by one of the surface functions.
260 vacuum: float
261 The thickness of the vacuum layer (in Angstrom).
262 """
263 uc = atoms.get_cell()
264 normal = np.cross(uc[0], uc[1])
265 costheta = np.dot(normal, uc[2]) / np.sqrt(np.dot(normal, normal) *
266 np.dot(uc[2], uc[2]))
267 length = np.sqrt(np.dot(uc[2], uc[2]))
268 newlength = length + vacuum / costheta
269 uc[2] *= newlength / length
270 atoms.set_cell(uc)
273def create_tags(size) -> np.ndarray:
274 """ Function to create layer tags. """
275 # tag atoms by layer
276 # create blocks of descending integers of length size[0]*size[1]
277 return np.arange(size[2], 0, -1).repeat(size[0] * size[1])
280def _surface(symbol, structure, face, size, a, c, vacuum, periodic,
281 orthogonal=True):
282 """Function to build often used surfaces.
284 Don't call this function directly - use fcc100, fcc110, bcc111, ..."""
286 Z = atomic_numbers[symbol]
288 if a is None:
289 sym = reference_states[Z]['symmetry']
290 if sym != structure:
291 raise ValueError(
292 f"Can't guess lattice constant for {structure}-{symbol}!")
293 a = reference_states[Z]['a']
295 if structure == 'hcp' and c is None:
296 if reference_states[Z]['symmetry'] == 'hcp':
297 c = reference_states[Z]['c/a'] * a
298 else:
299 c = sqrt(8 / 3.0) * a
301 positions = np.empty((size[2], size[1], size[0], 3))
302 positions[..., 0] = np.arange(size[0]).reshape((1, 1, -1))
303 positions[..., 1] = np.arange(size[1]).reshape((1, -1, 1))
304 positions[..., 2] = np.arange(size[2]).reshape((-1, 1, 1))
306 numbers = np.ones(size[0] * size[1] * size[2], int) * Z
308 slab = Atoms(numbers,
309 tags=create_tags(size),
310 pbc=(True, True, periodic),
311 cell=size)
313 surface_cell = None
314 sites = {'ontop': (0, 0)}
315 surf = structure + face
316 if surf == 'fcc100':
317 cell = (sqrt(0.5), sqrt(0.5), 0.5)
318 positions[-2::-2, ..., :2] += 0.5
319 sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)})
320 elif surf == 'diamond100':
321 cell = (sqrt(0.5), sqrt(0.5), 0.5 / 2)
322 positions[-4::-4, ..., :2] += (0.5, 0.5)
323 positions[-3::-4, ..., :2] += (0.0, 0.5)
324 positions[-2::-4, ..., :2] += (0.0, 0.0)
325 positions[-1::-4, ..., :2] += (0.5, 0.0)
326 elif surf == 'fcc110':
327 cell = (1.0, sqrt(0.5), sqrt(0.125))
328 positions[-2::-2, ..., :2] += 0.5
329 sites.update({'hollow': (0.5, 0.5), 'longbridge': (0.5, 0),
330 'shortbridge': (0, 0.5)})
331 elif surf == 'bcc100':
332 cell = (1.0, 1.0, 0.5)
333 positions[-2::-2, ..., :2] += 0.5
334 sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)})
335 else:
336 if orthogonal and size[1] % 2 == 1:
337 raise ValueError(("Can't make orthorhombic cell with size=%r. " %
338 (tuple(size),)) +
339 'Second number in size must be even.')
340 if surf == 'fcc111':
341 cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3))
342 if orthogonal:
343 positions[-1::-3, 1::2, :, 0] += 0.5
344 positions[-2::-3, 1::2, :, 0] += 0.5
345 positions[-3::-3, 1::2, :, 0] -= 0.5
346 positions[-2::-3, ..., :2] += (0.0, 2.0 / 3)
347 positions[-3::-3, ..., :2] += (0.5, 1.0 / 3)
348 else:
349 positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3)
350 positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3)
351 sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3),
352 'hcp': (2.0 / 3, 2.0 / 3)})
353 elif surf == 'diamond111':
354 cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3) / 2)
355 assert not orthogonal
356 positions[-1::-6, ..., :3] += (0.0, 0.0, 0.5)
357 positions[-2::-6, ..., :2] += (0.0, 0.0)
358 positions[-3::-6, ..., :3] += (-1.0 / 3, 2.0 / 3, 0.5)
359 positions[-4::-6, ..., :2] += (-1.0 / 3, 2.0 / 3)
360 positions[-5::-6, ..., :3] += (1.0 / 3, 1.0 / 3, 0.5)
361 positions[-6::-6, ..., :2] += (1.0 / 3, 1.0 / 3)
362 elif surf == 'hcp0001':
363 cell = (1.0, sqrt(0.75), 0.5 * c / a)
364 if orthogonal:
365 positions[:, 1::2, :, 0] += 0.5
366 positions[-2::-2, ..., :2] += (0.0, 2.0 / 3)
367 else:
368 positions[-2::-2, ..., :2] += (-1.0 / 3, 2.0 / 3)
369 sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3),
370 'hcp': (2.0 / 3, 2.0 / 3)})
371 elif surf == 'hcp10m10':
372 cell = (1.0, 0.5 * c / a, sqrt(0.75))
373 assert orthogonal
374 positions[-2::-2, ..., 0] += 0.5
375 positions[:, ::2, :, 2] += 2.0 / 3
376 elif surf == 'bcc110':
377 cell = (1.0, sqrt(0.5), sqrt(0.5))
378 if orthogonal:
379 positions[:, 1::2, :, 0] += 0.5
380 positions[-2::-2, ..., :2] += (0.0, 1.0)
381 else:
382 positions[-2::-2, ..., :2] += (-0.5, 1.0)
383 sites.update({'shortbridge': (0, 0.5),
384 'longbridge': (0.5, 0),
385 'hollow': (0.375, 0.25)})
386 elif surf == 'bcc111':
387 cell = (sqrt(2), sqrt(1.5), sqrt(3) / 6)
388 if orthogonal:
389 positions[-1::-3, 1::2, :, 0] += 0.5
390 positions[-2::-3, 1::2, :, 0] += 0.5
391 positions[-3::-3, 1::2, :, 0] -= 0.5
392 positions[-2::-3, ..., :2] += (0.0, 2.0 / 3)
393 positions[-3::-3, ..., :2] += (0.5, 1.0 / 3)
394 else:
395 positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3)
396 positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3)
397 sites.update({'hollow': (1.0 / 3, 1.0 / 3)})
398 else:
399 2 / 0
401 surface_cell = a * np.array([(cell[0], 0),
402 (cell[0] / 2, cell[1])])
403 if not orthogonal:
404 cell = np.array([(cell[0], 0, 0),
405 (cell[0] / 2, cell[1], 0),
406 (0, 0, cell[2])])
408 if surface_cell is None:
409 surface_cell = a * np.diag(cell[:2])
411 if isinstance(cell, tuple):
412 cell = np.diag(cell)
414 slab.set_positions(positions.reshape((-1, 3)))
415 slab.set_cell([a * v * n for v, n in zip(cell, size)], scale_atoms=True)
417 if not periodic:
418 slab.cell[2] = 0.0
420 if vacuum is not None:
421 slab.center(vacuum, axis=2)
423 if 'adsorbate_info' not in slab.info:
424 slab.info.update({'adsorbate_info': {}})
426 slab.info['adsorbate_info']['cell'] = surface_cell
427 slab.info['adsorbate_info']['sites'] = sites
428 return slab
431def fcc211(symbol, size, a=None, vacuum=None, orthogonal=True):
432 """FCC(211) surface.
434 Does not currently support special adsorption sites.
436 Currently only implemented for *orthogonal=True* with size specified
437 as (i, j, k), where i, j, and k are number of atoms in each direction.
438 i must be divisible by 3 to accommodate the step width.
439 """
440 if not orthogonal:
441 raise NotImplementedError('Only implemented for orthogonal '
442 'unit cells.')
443 if size[0] % 3 != 0:
444 raise NotImplementedError('First dimension of size must be '
445 'divisible by 3.')
446 atoms = FaceCenteredCubic(symbol,
447 directions=[[1, -1, -1],
448 [0, 2, -2],
449 [2, 1, 1]],
450 miller=(None, None, (2, 1, 1)),
451 latticeconstant=a,
452 size=(1, 1, 1),
453 pbc=True)
454 z = (size[2] + 1) // 2
455 atoms = atoms.repeat((size[0] // 3, size[1], z))
456 if size[2] % 2: # Odd: remove bottom layer and shrink cell.
457 remove_list = [atom.index for atom in atoms
458 if atom.z < atoms[1].z]
459 del atoms[remove_list]
460 dz = atoms[0].z
461 atoms.translate((0., 0., -dz))
462 atoms.cell[2][2] -= dz
464 atoms.cell[2] = 0.0
465 atoms.pbc[2] = False
466 if vacuum:
467 atoms.center(vacuum, axis=2)
469 # Renumber systematically from top down.
470 orders = [(atom.index, round(atom.x, 3), round(atom.y, 3),
471 -round(atom.z, 3), atom.index) for atom in atoms]
472 orders.sort(key=itemgetter(3, 1, 2))
473 newatoms = atoms.copy()
474 for index, order in enumerate(orders):
475 newatoms[index].position = atoms[order[0]].position.copy()
477 # Add empty 'sites' dictionary for consistency with other functions
478 newatoms.info['adsorbate_info'] = {'sites': {}}
479 return newatoms
482def mx2(formula='MoS2', kind='2H', a=3.18, thickness=3.19,
483 size=(1, 1, 1), vacuum=None):
484 """Create three-layer 2D materials with hexagonal structure.
486 This can be used for e.g. metal dichalcogenides :mol:`MX_2` 2D structures
487 such as :mol:`MoS_2`.
489 https://en.wikipedia.org/wiki/Transition_metal_dichalcogenide_monolayers
491 Parameters
492 ----------
493 kind : {'2H', '1T'}, default: '2H'
495 - '2H': mirror-plane symmetry
496 - '1T': inversion symmetry
497 """
498 if kind == '2H':
499 basis = [(0, 0, 0),
500 (2 / 3, 1 / 3, 0.5 * thickness),
501 (2 / 3, 1 / 3, -0.5 * thickness)]
502 elif kind == '1T':
503 basis = [(0, 0, 0),
504 (2 / 3, 1 / 3, 0.5 * thickness),
505 (1 / 3, 2 / 3, -0.5 * thickness)]
506 else:
507 raise ValueError('Structure not recognized:', kind)
509 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
511 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
512 atoms.set_scaled_positions(basis)
513 if vacuum is not None:
514 atoms.center(vacuum, axis=2)
515 atoms = atoms.repeat(size)
516 return atoms
519def graphene(formula='C2', a=2.460, thickness=0.0,
520 size=(1, 1, 1), vacuum=None):
521 """Create a graphene monolayer structure.
523 Parameters
524 ----------
525 thickness : float, default: 0.0
526 Thickness of the layer; maybe for a buckled structure like silicene.
527 """
528 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
529 basis = [[0, 0, -0.5 * thickness], [2 / 3, 1 / 3, 0.5 * thickness]]
530 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
531 atoms.set_scaled_positions(basis)
532 if vacuum is not None:
533 atoms.center(vacuum, axis=2)
534 atoms = atoms.repeat(size)
535 return atoms
538def _all_surface_functions():
539 # Convenient for debugging.
540 d = {
541 func.__name__: func
542 for func in [
543 fcc100,
544 fcc110,
545 bcc100,
546 bcc110,
547 bcc111,
548 fcc111,
549 hcp0001,
550 hcp10m10,
551 diamond100,
552 diamond111,
553 fcc111,
554 mx2,
555 graphene,
556 ]
557 }
558 return d