Coverage for /builds/kinetik161/ase/ase/build/surface.py: 84.23%
241 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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 reference_states, atomic_numbers
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 _surface(symbol, structure, face, size, a, c, vacuum, periodic,
272 orthogonal=True):
273 """Function to build often used surfaces.
275 Don't call this function directly - use fcc100, fcc110, bcc111, ..."""
277 Z = atomic_numbers[symbol]
279 if a is None:
280 sym = reference_states[Z]['symmetry']
281 if sym != structure:
282 raise ValueError(
283 f"Can't guess lattice constant for {structure}-{symbol}!")
284 a = reference_states[Z]['a']
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
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))
297 numbers = np.ones(size[0] * size[1] * size[2], int) * Z
299 tags = np.empty((size[2], size[1], size[0]), int)
300 tags[:] = np.arange(size[2], 0, -1).reshape((-1, 1, 1))
302 slab = Atoms(numbers,
303 tags=tags.ravel(),
304 pbc=(True, True, periodic),
305 cell=size)
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
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])])
402 if surface_cell is None:
403 surface_cell = a * np.diag(cell[:2])
405 if isinstance(cell, tuple):
406 cell = np.diag(cell)
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)
411 if not periodic:
412 slab.cell[2] = 0.0
414 if vacuum is not None:
415 slab.center(vacuum, axis=2)
417 if 'adsorbate_info' not in slab.info:
418 slab.info.update({'adsorbate_info': {}})
420 slab.info['adsorbate_info']['cell'] = surface_cell
421 slab.info['adsorbate_info']['sites'] = sites
422 return slab
425def fcc211(symbol, size, a=None, vacuum=None, orthogonal=True):
426 """FCC(211) surface.
428 Does not currently support special adsorption sites.
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
458 atoms.cell[2] = 0.0
459 atoms.pbc[2] = False
460 if vacuum:
461 atoms.center(vacuum, axis=2)
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()
471 # Add empty 'sites' dictionary for consistency with other functions
472 newatoms.info['adsorbate_info'] = {'sites': {}}
473 return newatoms
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.
480 This can be used for e.g. metal dichalcogenides :mol:`MX_2` 2D structures
481 such as :mol:`MoS_2`.
483 https://en.wikipedia.org/wiki/Transition_metal_dichalcogenide_monolayers
485 Parameters
486 ----------
487 kind : {'2H', '1T'}, default: '2H'
489 - '2H': mirror-plane symmetry
490 - '1T': inversion symmetry
491 """
492 if kind == '2H':
493 basis = [(0, 0, 0),
494 (2 / 3, 1 / 3, 0.5 * thickness),
495 (2 / 3, 1 / 3, -0.5 * thickness)]
496 elif kind == '1T':
497 basis = [(0, 0, 0),
498 (2 / 3, 1 / 3, 0.5 * thickness),
499 (1 / 3, 2 / 3, -0.5 * thickness)]
500 else:
501 raise ValueError('Structure not recognized:', kind)
503 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
505 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
506 atoms.set_scaled_positions(basis)
507 if vacuum is not None:
508 atoms.center(vacuum, axis=2)
509 atoms = atoms.repeat(size)
510 return atoms
513def graphene(formula='C2', a=2.460, thickness=0.0,
514 size=(1, 1, 1), vacuum=None):
515 """Create a graphene monolayer structure.
517 Parameters
518 ----------
519 thickness : float, default: 0.0
520 Thickness of the layer; maybe for a buckled structure like silicene.
521 """
522 cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
523 basis = [[0, 0, -0.5 * thickness], [2 / 3, 1 / 3, 0.5 * thickness]]
524 atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
525 atoms.set_scaled_positions(basis)
526 if vacuum is not None:
527 atoms.center(vacuum, axis=2)
528 atoms = atoms.repeat(size)
529 return atoms
532def _all_surface_functions():
533 # Convenient for debugging.
534 d = {}
535 for func in [fcc100, fcc110, bcc100, bcc110, bcc111, fcc111, hcp0001,
536 hcp10m10, diamond100, diamond111, fcc111, mx2, graphene]:
537 d[func.__name__] = func
538 return d