Coverage for /builds/debichem-team/python-ase/ase/build/surface.py: 84.10%
239 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1"""Helper functions for creating the most common surfaces and related tasks.
3The helper functions can create the most common low-index surfaces,
4add vacuum layers and add adsorbates.
6"""
8from math import sqrt
9from operator import itemgetter
11import numpy as np
13from ase.atom import Atom
14from ase.atoms import Atoms
15from ase.data import atomic_numbers, reference_states
16from ase.lattice.cubic import FaceCenteredCubic
19def fcc100(symbol, size, a=None, vacuum=None, orthogonal=True,
20 periodic=False):
21 """FCC(100) surface.
23 Supported special adsorption sites: 'ontop', 'bridge', 'hollow'."""
24 if not orthogonal:
25 raise NotImplementedError("Can't do non-orthogonal cell yet!")
27 return _surface(symbol, 'fcc', '100', size, a, None, vacuum,
28 periodic=periodic,
29 orthogonal=orthogonal)
32def fcc110(symbol, size, a=None, vacuum=None, orthogonal=True,
33 periodic=False):
34 """FCC(110) surface.
36 Supported special adsorption sites: 'ontop', 'longbridge',
37 'shortbridge', 'hollow'."""
38 if not orthogonal:
39 raise NotImplementedError("Can't do non-orthogonal cell yet!")
41 return _surface(symbol, 'fcc', '110', size, a, None, vacuum,
42 periodic=periodic,
43 orthogonal=orthogonal)
46def bcc100(symbol, size, a=None, vacuum=None, orthogonal=True,
47 periodic=False):
48 """BCC(100) surface.
50 Supported special adsorption sites: 'ontop', 'bridge', 'hollow'."""
51 if not orthogonal:
52 raise NotImplementedError("Can't do non-orthogonal cell yet!")
54 return _surface(symbol, 'bcc', '100', size, a, None, vacuum,
55 periodic=periodic,
56 orthogonal=orthogonal)
59def bcc110(symbol, size, a=None, vacuum=None, orthogonal=False,
60 periodic=False):
61 """BCC(110) surface.
63 Supported special adsorption sites: 'ontop', 'longbridge',
64 'shortbridge', 'hollow'.
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)
73def bcc111(symbol, size, a=None, vacuum=None, orthogonal=False,
74 periodic=False):
75 """BCC(111) surface.
77 Supported special adsorption sites: 'ontop'.
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)
86def fcc111(symbol, size, a=None, vacuum=None, orthogonal=False,
87 periodic=False):
88 """FCC(111) surface.
90 Supported special adsorption sites: 'ontop', 'bridge', 'fcc' and 'hcp'.
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)
99def hcp0001(symbol, size, a=None, c=None, vacuum=None, orthogonal=False,
100 periodic=False):
101 """HCP(0001) surface.
103 Supported special adsorption sites: 'ontop', 'bridge', 'fcc' and 'hcp'.
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)
112def hcp10m10(symbol, size, a=None, c=None, vacuum=None, orthogonal=True,
113 periodic=False):
114 """HCP(10m10) surface.
116 Supported special adsorption sites: 'ontop'.
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!")
122 return _surface(symbol, 'hcp', '10m10', size, a, c, vacuum,
123 periodic=periodic,
124 orthogonal=orthogonal)
127def diamond100(symbol, size, a=None, vacuum=None, orthogonal=True,
128 periodic=False):
129 """DIAMOND(100) surface.
131 Supported special adsorption sites: 'ontop'."""
132 if not orthogonal:
133 raise NotImplementedError("Can't do non-orthogonal cell yet!")
135 return _surface(symbol, 'diamond', '100', size, a, None, vacuum,
136 periodic=periodic,
137 orthogonal=orthogonal)
140def diamond111(symbol, size, a=None, vacuum=None, orthogonal=False,
141 periodic=False):
142 """DIAMOND(111) surface.
144 Supported special adsorption sites: 'ontop'."""
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)
153def add_adsorbate(slab, adsorbate, height, position=(0, 0), offset=None,
154 mol_index=0):
155 """Add an adsorbate to a surface.
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).
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.
168 This function can be called multiple times to add more than one
169 adsorbate.
171 Parameters:
173 slab: The surface onto which the adsorbate should be added.
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).
180 height: Height above the surface.
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).
186 offset (default: None): Offsets the adsorbate by a number of unit
187 cells. Mostly useful when adding more than one adsorbate.
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.
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.
198 """
199 info = slab.info.get('adsorbate_info', {})
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)
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(f'Adsorption site {position} not supported.')
214 spos += info['sites'][position]
215 else:
216 pos += position
218 if 'cell' in info:
219 cell = info['cell']
220 else:
221 cell = slab.get_cell()[:2, :2]
223 pos += np.dot(spos, cell)
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)])
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
244 # Move adsorbate into position
245 ads.translate([pos[0], pos[1], z] - ads.positions[mol_index])
247 # Attach the adsorbate
248 slab.extend(ads)
251def add_vacuum(atoms, vacuum):
252 """Add vacuum layer to the atoms.
254 Parameters:
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)
271def create_tags(size) -> np.ndarray:
272 """ Function to create layer tags. """
273 # tag atoms by layer
274 # create blocks of descending integers of length size[0]*size[1]
275 return np.arange(size[2], 0, -1).repeat(size[0] * size[1])
278def _surface(symbol, structure, face, size, a, c, vacuum, periodic,
279 orthogonal=True):
280 """Function to build often used surfaces.
282 Don't call this function directly - use fcc100, fcc110, bcc111, ..."""
284 Z = atomic_numbers[symbol]
286 if a is None:
287 sym = reference_states[Z]['symmetry']
288 if sym != structure:
289 raise ValueError(
290 f"Can't guess lattice constant for {structure}-{symbol}!")
291 a = reference_states[Z]['a']
293 if structure == 'hcp' and c is None:
294 if reference_states[Z]['symmetry'] == 'hcp':
295 c = reference_states[Z]['c/a'] * a
296 else:
297 c = sqrt(8 / 3.0) * a
299 positions = np.empty((size[2], size[1], size[0], 3))
300 positions[..., 0] = np.arange(size[0]).reshape((1, 1, -1))
301 positions[..., 1] = np.arange(size[1]).reshape((1, -1, 1))
302 positions[..., 2] = np.arange(size[2]).reshape((-1, 1, 1))
304 numbers = np.ones(size[0] * size[1] * size[2], int) * Z
306 slab = Atoms(numbers,
307 tags=create_tags(size),
308 pbc=(True, True, periodic),
309 cell=size)
311 surface_cell = None
312 sites = {'ontop': (0, 0)}
313 surf = structure + face
314 if surf == 'fcc100':
315 cell = (sqrt(0.5), sqrt(0.5), 0.5)
316 positions[-2::-2, ..., :2] += 0.5
317 sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)})
318 elif surf == 'diamond100':
319 cell = (sqrt(0.5), sqrt(0.5), 0.5 / 2)
320 positions[-4::-4, ..., :2] += (0.5, 0.5)
321 positions[-3::-4, ..., :2] += (0.0, 0.5)
322 positions[-2::-4, ..., :2] += (0.0, 0.0)
323 positions[-1::-4, ..., :2] += (0.5, 0.0)
324 elif surf == 'fcc110':
325 cell = (1.0, sqrt(0.5), sqrt(0.125))
326 positions[-2::-2, ..., :2] += 0.5
327 sites.update({'hollow': (0.5, 0.5), 'longbridge': (0.5, 0),
328 'shortbridge': (0, 0.5)})
329 elif surf == 'bcc100':
330 cell = (1.0, 1.0, 0.5)
331 positions[-2::-2, ..., :2] += 0.5
332 sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)})
333 else:
334 if orthogonal and size[1] % 2 == 1:
335 raise ValueError(("Can't make orthorhombic cell with size=%r. " %
336 (tuple(size),)) +
337 'Second number in size must be even.')
338 if surf == 'fcc111':
339 cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3))
340 if orthogonal:
341 positions[-1::-3, 1::2, :, 0] += 0.5
342 positions[-2::-3, 1::2, :, 0] += 0.5
343 positions[-3::-3, 1::2, :, 0] -= 0.5
344 positions[-2::-3, ..., :2] += (0.0, 2.0 / 3)
345 positions[-3::-3, ..., :2] += (0.5, 1.0 / 3)
346 else:
347 positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3)
348 positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3)
349 sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3),
350 'hcp': (2.0 / 3, 2.0 / 3)})
351 elif surf == 'diamond111':
352 cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3) / 2)
353 assert not orthogonal
354 positions[-1::-6, ..., :3] += (0.0, 0.0, 0.5)
355 positions[-2::-6, ..., :2] += (0.0, 0.0)
356 positions[-3::-6, ..., :3] += (-1.0 / 3, 2.0 / 3, 0.5)
357 positions[-4::-6, ..., :2] += (-1.0 / 3, 2.0 / 3)
358 positions[-5::-6, ..., :3] += (1.0 / 3, 1.0 / 3, 0.5)
359 positions[-6::-6, ..., :2] += (1.0 / 3, 1.0 / 3)
360 elif surf == 'hcp0001':
361 cell = (1.0, sqrt(0.75), 0.5 * c / a)
362 if orthogonal:
363 positions[:, 1::2, :, 0] += 0.5
364 positions[-2::-2, ..., :2] += (0.0, 2.0 / 3)
365 else:
366 positions[-2::-2, ..., :2] += (-1.0 / 3, 2.0 / 3)
367 sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3),
368 'hcp': (2.0 / 3, 2.0 / 3)})
369 elif surf == 'hcp10m10':
370 cell = (1.0, 0.5 * c / a, sqrt(0.75))
371 assert orthogonal
372 positions[-2::-2, ..., 0] += 0.5
373 positions[:, ::2, :, 2] += 2.0 / 3
374 elif surf == 'bcc110':
375 cell = (1.0, sqrt(0.5), sqrt(0.5))
376 if orthogonal:
377 positions[:, 1::2, :, 0] += 0.5
378 positions[-2::-2, ..., :2] += (0.0, 1.0)
379 else:
380 positions[-2::-2, ..., :2] += (-0.5, 1.0)
381 sites.update({'shortbridge': (0, 0.5),
382 'longbridge': (0.5, 0),
383 'hollow': (0.375, 0.25)})
384 elif surf == 'bcc111':
385 cell = (sqrt(2), sqrt(1.5), sqrt(3) / 6)
386 if orthogonal:
387 positions[-1::-3, 1::2, :, 0] += 0.5
388 positions[-2::-3, 1::2, :, 0] += 0.5
389 positions[-3::-3, 1::2, :, 0] -= 0.5
390 positions[-2::-3, ..., :2] += (0.0, 2.0 / 3)
391 positions[-3::-3, ..., :2] += (0.5, 1.0 / 3)
392 else:
393 positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3)
394 positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3)
395 sites.update({'hollow': (1.0 / 3, 1.0 / 3)})
396 else:
397 2 / 0
399 surface_cell = a * np.array([(cell[0], 0),
400 (cell[0] / 2, cell[1])])
401 if not orthogonal:
402 cell = np.array([(cell[0], 0, 0),
403 (cell[0] / 2, cell[1], 0),
404 (0, 0, cell[2])])
406 if surface_cell is None:
407 surface_cell = a * np.diag(cell[:2])
409 if isinstance(cell, tuple):
410 cell = np.diag(cell)
412 slab.set_positions(positions.reshape((-1, 3)))
413 slab.set_cell([a * v * n for v, n in zip(cell, size)], scale_atoms=True)
415 if not periodic:
416 slab.cell[2] = 0.0
418 if vacuum is not None:
419 slab.center(vacuum, axis=2)
421 if 'adsorbate_info' not in slab.info:
422 slab.info.update({'adsorbate_info': {}})
424 slab.info['adsorbate_info']['cell'] = surface_cell
425 slab.info['adsorbate_info']['sites'] = sites
426 return slab
429def fcc211(symbol, size, a=None, vacuum=None, orthogonal=True):
430 """FCC(211) surface.
432 Does not currently support special adsorption sites.
434 Currently only implemented for *orthogonal=True* with size specified
435 as (i, j, k), where i, j, and k are number of atoms in each direction.
436 i must be divisible by 3 to accommodate the step width.
437 """
438 if not orthogonal:
439 raise NotImplementedError('Only implemented for orthogonal '
440 'unit cells.')
441 if size[0] % 3 != 0:
442 raise NotImplementedError('First dimension of size must be '
443 'divisible by 3.')
444 atoms = FaceCenteredCubic(symbol,
445 directions=[[1, -1, -1],
446 [0, 2, -2],
447 [2, 1, 1]],
448 miller=(None, None, (2, 1, 1)),
449 latticeconstant=a,
450 size=(1, 1, 1),
451 pbc=True)
452 z = (size[2] + 1) // 2
453 atoms = atoms.repeat((size[0] // 3, size[1], z))
454 if size[2] % 2: # Odd: remove bottom layer and shrink cell.
455 remove_list = [atom.index for atom in atoms
456 if atom.z < atoms[1].z]
457 del atoms[remove_list]
458 dz = atoms[0].z
459 atoms.translate((0., 0., -dz))
460 atoms.cell[2][2] -= dz
462 atoms.cell[2] = 0.0
463 atoms.pbc[2] = False
464 if vacuum:
465 atoms.center(vacuum, axis=2)
467 # Renumber systematically from top down.
468 orders = [(atom.index, round(atom.x, 3), round(atom.y, 3),
469 -round(atom.z, 3), atom.index) for atom in atoms]
470 orders.sort(key=itemgetter(3, 1, 2))
471 newatoms = atoms.copy()
472 for index, order in enumerate(orders):
473 newatoms[index].position = atoms[order[0]].position.copy()
475 # Add empty 'sites' dictionary for consistency with other functions
476 newatoms.info['adsorbate_info'] = {'sites': {}}
477 return newatoms
480def mx2(formula='MoS2', kind='2H', a=3.18, thickness=3.19,
481 size=(1, 1, 1), vacuum=None):
482 """Create three-layer 2D materials with hexagonal structure.
484 This can be used for e.g. metal dichalcogenides :mol:`MX_2` 2D structures
485 such as :mol:`MoS_2`.
487 https://en.wikipedia.org/wiki/Transition_metal_dichalcogenide_monolayers
489 Parameters
490 ----------
491 kind : {'2H', '1T'}, default: '2H'
493 - '2H': mirror-plane symmetry
494 - '1T': inversion symmetry
495 """
496 if kind == '2H':
497 basis = [(0, 0, 0),
498 (2 / 3, 1 / 3, 0.5 * thickness),
499 (2 / 3, 1 / 3, -0.5 * thickness)]
500 elif kind == '1T':
501 basis = [(0, 0, 0),
502 (2 / 3, 1 / 3, 0.5 * thickness),
503 (1 / 3, 2 / 3, -0.5 * thickness)]
504 else:
505 raise ValueError('Structure not recognized:', kind)
507 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
509 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
510 atoms.set_scaled_positions(basis)
511 if vacuum is not None:
512 atoms.center(vacuum, axis=2)
513 atoms = atoms.repeat(size)
514 return atoms
517def graphene(formula='C2', a=2.460, thickness=0.0,
518 size=(1, 1, 1), vacuum=None):
519 """Create a graphene monolayer structure.
521 Parameters
522 ----------
523 thickness : float, default: 0.0
524 Thickness of the layer; maybe for a buckled structure like silicene.
525 """
526 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
527 basis = [[0, 0, -0.5 * thickness], [2 / 3, 1 / 3, 0.5 * thickness]]
528 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
529 atoms.set_scaled_positions(basis)
530 if vacuum is not None:
531 atoms.center(vacuum, axis=2)
532 atoms = atoms.repeat(size)
533 return atoms
536def _all_surface_functions():
537 # Convenient for debugging.
538 d = {
539 func.__name__: func
540 for func in [
541 fcc100,
542 fcc110,
543 bcc100,
544 bcc110,
545 bcc111,
546 fcc111,
547 hcp0001,
548 hcp10m10,
549 diamond100,
550 diamond111,
551 fcc111,
552 mx2,
553 graphene,
554 ]
555 }
556 return d