"""Helper functions for creating the most common surfaces and related tasks.
The helper functions can create the most common low-index surfaces,
add vacuum layers and add adsorbates.
"""
from math import sqrt
from operator import itemgetter
import numpy as np
from ase.atom import Atom
from ase.atoms import Atoms
from ase.data import reference_states, atomic_numbers
from ase.lattice.cubic import FaceCenteredCubic
[docs]def fcc100(symbol, size, a=None, vacuum=None, orthogonal=True,
periodic=False):
"""FCC(100) surface.
Supported special adsorption sites: 'ontop', 'bridge', 'hollow'."""
if not orthogonal:
raise NotImplementedError("Can't do non-orthogonal cell yet!")
return _surface(symbol, 'fcc', '100', size, a, None, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def fcc110(symbol, size, a=None, vacuum=None, orthogonal=True,
periodic=False):
"""FCC(110) surface.
Supported special adsorption sites: 'ontop', 'longbridge',
'shortbridge', 'hollow'."""
if not orthogonal:
raise NotImplementedError("Can't do non-orthogonal cell yet!")
return _surface(symbol, 'fcc', '110', size, a, None, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def bcc100(symbol, size, a=None, vacuum=None, orthogonal=True,
periodic=False):
"""BCC(100) surface.
Supported special adsorption sites: 'ontop', 'bridge', 'hollow'."""
if not orthogonal:
raise NotImplementedError("Can't do non-orthogonal cell yet!")
return _surface(symbol, 'bcc', '100', size, a, None, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def bcc110(symbol, size, a=None, vacuum=None, orthogonal=False,
periodic=False):
"""BCC(110) surface.
Supported special adsorption sites: 'ontop', 'longbridge',
'shortbridge', 'hollow'.
Use *orthogonal=True* to get an orthogonal unit cell - works only
for size=(i,j,k) with j even."""
return _surface(symbol, 'bcc', '110', size, a, None, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def bcc111(symbol, size, a=None, vacuum=None, orthogonal=False,
periodic=False):
"""BCC(111) surface.
Supported special adsorption sites: 'ontop'.
Use *orthogonal=True* to get an orthogonal unit cell - works only
for size=(i,j,k) with j even."""
return _surface(symbol, 'bcc', '111', size, a, None, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def fcc111(symbol, size, a=None, vacuum=None, orthogonal=False,
periodic=False):
"""FCC(111) surface.
Supported special adsorption sites: 'ontop', 'bridge', 'fcc' and 'hcp'.
Use *orthogonal=True* to get an orthogonal unit cell - works only
for size=(i,j,k) with j even."""
return _surface(symbol, 'fcc', '111', size, a, None, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def hcp0001(symbol, size, a=None, c=None, vacuum=None, orthogonal=False,
periodic=False):
"""HCP(0001) surface.
Supported special adsorption sites: 'ontop', 'bridge', 'fcc' and 'hcp'.
Use *orthogonal=True* to get an orthogonal unit cell - works only
for size=(i,j,k) with j even."""
return _surface(symbol, 'hcp', '0001', size, a, c, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def hcp10m10(symbol, size, a=None, c=None, vacuum=None, orthogonal=True,
periodic=False):
"""HCP(10m10) surface.
Supported special adsorption sites: 'ontop'.
Works only for size=(i,j,k) with j even."""
if not orthogonal:
raise NotImplementedError("Can't do non-orthogonal cell yet!")
return _surface(symbol, 'hcp', '10m10', size, a, c, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def diamond100(symbol, size, a=None, vacuum=None, orthogonal=True,
periodic=False):
"""DIAMOND(100) surface.
Supported special adsorption sites: 'ontop'."""
if not orthogonal:
raise NotImplementedError("Can't do non-orthogonal cell yet!")
return _surface(symbol, 'diamond', '100', size, a, None, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def diamond111(symbol, size, a=None, vacuum=None, orthogonal=False,
periodic=False):
"""DIAMOND(111) surface.
Supported special adsorption sites: 'ontop'."""
if orthogonal:
raise NotImplementedError("Can't do orthogonal cell yet!")
return _surface(symbol, 'diamond', '111', size, a, None, vacuum,
periodic=periodic,
orthogonal=orthogonal)
[docs]def add_adsorbate(slab, adsorbate, height, position=(0, 0), offset=None,
mol_index=0):
"""Add an adsorbate to a surface.
This function adds an adsorbate to a slab. If the slab is
produced by one of the utility functions in ase.build, it
is possible to specify the position of the adsorbate by a keyword
(the supported keywords depend on which function was used to
create the slab).
If the adsorbate is a molecule, the atom indexed by the mol_index
optional argument is positioned on top of the adsorption position
on the surface, and it is the responsibility of the user to orient
the adsorbate in a sensible way.
This function can be called multiple times to add more than one
adsorbate.
Parameters:
slab: The surface onto which the adsorbate should be added.
adsorbate: The adsorbate. Must be one of the following three types:
A string containing the chemical symbol for a single atom.
An atom object.
An atoms object (for a molecular adsorbate).
height: Height above the surface.
position: The x-y position of the adsorbate, either as a tuple of
two numbers or as a keyword (if the surface is produced by one
of the functions in ase.build).
offset (default: None): Offsets the adsorbate by a number of unit
cells. Mostly useful when adding more than one adsorbate.
mol_index (default: 0): If the adsorbate is a molecule, index of
the atom to be positioned above the location specified by the
position argument.
Note *position* is given in absolute xy coordinates (or as
a keyword), whereas offset is specified in unit cells. This
can be used to give the positions in units of the unit cell by
using *offset* instead.
"""
info = slab.info.get('adsorbate_info', {})
pos = np.array([0.0, 0.0]) # (x, y) part
spos = np.array([0.0, 0.0]) # part relative to unit cell
if offset is not None:
spos += np.asarray(offset, float)
if isinstance(position, str):
# A site-name:
if 'sites' not in info:
raise TypeError('If the atoms are not made by an ' +
'ase.build function, ' +
'position cannot be a name.')
if position not in info['sites']:
raise TypeError('Adsorption site %s not supported.' % position)
spos += info['sites'][position]
else:
pos += position
if 'cell' in info:
cell = info['cell']
else:
cell = slab.get_cell()[:2, :2]
pos += np.dot(spos, cell)
# Convert the adsorbate to an Atoms object
if isinstance(adsorbate, Atoms):
ads = adsorbate
elif isinstance(adsorbate, Atom):
ads = Atoms([adsorbate])
else:
# Assume it is a string representing a single Atom
ads = Atoms([Atom(adsorbate)])
# Get the z-coordinate:
if 'top layer atom index' in info:
a = info['top layer atom index']
else:
a = slab.positions[:, 2].argmax()
if 'adsorbate_info' not in slab.info:
slab.info['adsorbate_info'] = {}
slab.info['adsorbate_info']['top layer atom index'] = a
z = slab.positions[a, 2] + height
# Move adsorbate into position
ads.translate([pos[0], pos[1], z] - ads.positions[mol_index])
# Attach the adsorbate
slab.extend(ads)
[docs]def add_vacuum(atoms, vacuum):
"""Add vacuum layer to the atoms.
Parameters:
atoms: Atoms object
Most likely created by one of the surface functions.
vacuum: float
The thickness of the vacuum layer (in Angstrom).
"""
uc = atoms.get_cell()
normal = np.cross(uc[0], uc[1])
costheta = np.dot(normal, uc[2]) / np.sqrt(np.dot(normal, normal) *
np.dot(uc[2], uc[2]))
length = np.sqrt(np.dot(uc[2], uc[2]))
newlength = length + vacuum / costheta
uc[2] *= newlength / length
atoms.set_cell(uc)
def _surface(symbol, structure, face, size, a, c, vacuum, periodic,
orthogonal=True):
"""Function to build often used surfaces.
Don't call this function directly - use fcc100, fcc110, bcc111, ..."""
Z = atomic_numbers[symbol]
if a is None:
sym = reference_states[Z]['symmetry']
if sym != structure:
raise ValueError("Can't guess lattice constant for %s-%s!" %
(structure, symbol))
a = reference_states[Z]['a']
if structure == 'hcp' and c is None:
if reference_states[Z]['symmetry'] == 'hcp':
c = reference_states[Z]['c/a'] * a
else:
c = sqrt(8 / 3.0) * a
positions = np.empty((size[2], size[1], size[0], 3))
positions[..., 0] = np.arange(size[0]).reshape((1, 1, -1))
positions[..., 1] = np.arange(size[1]).reshape((1, -1, 1))
positions[..., 2] = np.arange(size[2]).reshape((-1, 1, 1))
numbers = np.ones(size[0] * size[1] * size[2], int) * Z
tags = np.empty((size[2], size[1], size[0]), int)
tags[:] = np.arange(size[2], 0, -1).reshape((-1, 1, 1))
slab = Atoms(numbers,
tags=tags.ravel(),
pbc=(True, True, periodic),
cell=size)
surface_cell = None
sites = {'ontop': (0, 0)}
surf = structure + face
if surf == 'fcc100':
cell = (sqrt(0.5), sqrt(0.5), 0.5)
positions[-2::-2, ..., :2] += 0.5
sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)})
elif surf == 'diamond100':
cell = (sqrt(0.5), sqrt(0.5), 0.5 / 2)
positions[-4::-4, ..., :2] += (0.5, 0.5)
positions[-3::-4, ..., :2] += (0.0, 0.5)
positions[-2::-4, ..., :2] += (0.0, 0.0)
positions[-1::-4, ..., :2] += (0.5, 0.0)
elif surf == 'fcc110':
cell = (1.0, sqrt(0.5), sqrt(0.125))
positions[-2::-2, ..., :2] += 0.5
sites.update({'hollow': (0.5, 0.5), 'longbridge': (0.5, 0),
'shortbridge': (0, 0.5)})
elif surf == 'bcc100':
cell = (1.0, 1.0, 0.5)
positions[-2::-2, ..., :2] += 0.5
sites.update({'hollow': (0.5, 0.5), 'bridge': (0.5, 0)})
else:
if orthogonal and size[1] % 2 == 1:
raise ValueError(("Can't make orthorhombic cell with size=%r. " %
(tuple(size),)) +
'Second number in size must be even.')
if surf == 'fcc111':
cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3))
if orthogonal:
positions[-1::-3, 1::2, :, 0] += 0.5
positions[-2::-3, 1::2, :, 0] += 0.5
positions[-3::-3, 1::2, :, 0] -= 0.5
positions[-2::-3, ..., :2] += (0.0, 2.0 / 3)
positions[-3::-3, ..., :2] += (0.5, 1.0 / 3)
else:
positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3)
positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3)
sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3),
'hcp': (2.0 / 3, 2.0 / 3)})
elif surf == 'diamond111':
cell = (sqrt(0.5), sqrt(0.375), 1 / sqrt(3) / 2)
assert not orthogonal
positions[-1::-6, ..., :3] += (0.0, 0.0, 0.5)
positions[-2::-6, ..., :2] += (0.0, 0.0)
positions[-3::-6, ..., :3] += (-1.0 / 3, 2.0 / 3, 0.5)
positions[-4::-6, ..., :2] += (-1.0 / 3, 2.0 / 3)
positions[-5::-6, ..., :3] += (1.0 / 3, 1.0 / 3, 0.5)
positions[-6::-6, ..., :2] += (1.0 / 3, 1.0 / 3)
elif surf == 'hcp0001':
cell = (1.0, sqrt(0.75), 0.5 * c / a)
if orthogonal:
positions[:, 1::2, :, 0] += 0.5
positions[-2::-2, ..., :2] += (0.0, 2.0 / 3)
else:
positions[-2::-2, ..., :2] += (-1.0 / 3, 2.0 / 3)
sites.update({'bridge': (0.5, 0), 'fcc': (1.0 / 3, 1.0 / 3),
'hcp': (2.0 / 3, 2.0 / 3)})
elif surf == 'hcp10m10':
cell = (1.0, 0.5 * c / a, sqrt(0.75))
assert orthogonal
positions[-2::-2, ..., 0] += 0.5
positions[:, ::2, :, 2] += 2.0 / 3
elif surf == 'bcc110':
cell = (1.0, sqrt(0.5), sqrt(0.5))
if orthogonal:
positions[:, 1::2, :, 0] += 0.5
positions[-2::-2, ..., :2] += (0.0, 1.0)
else:
positions[-2::-2, ..., :2] += (-0.5, 1.0)
sites.update({'shortbridge': (0, 0.5),
'longbridge': (0.5, 0),
'hollow': (0.375, 0.25)})
elif surf == 'bcc111':
cell = (sqrt(2), sqrt(1.5), sqrt(3) / 6)
if orthogonal:
positions[-1::-3, 1::2, :, 0] += 0.5
positions[-2::-3, 1::2, :, 0] += 0.5
positions[-3::-3, 1::2, :, 0] -= 0.5
positions[-2::-3, ..., :2] += (0.0, 2.0 / 3)
positions[-3::-3, ..., :2] += (0.5, 1.0 / 3)
else:
positions[-2::-3, ..., :2] += (-1.0 / 3, 2.0 / 3)
positions[-3::-3, ..., :2] += (1.0 / 3, 1.0 / 3)
sites.update({'hollow': (1.0 / 3, 1.0 / 3)})
else:
2 / 0
surface_cell = a * np.array([(cell[0], 0),
(cell[0] / 2, cell[1])])
if not orthogonal:
cell = np.array([(cell[0], 0, 0),
(cell[0] / 2, cell[1], 0),
(0, 0, cell[2])])
if surface_cell is None:
surface_cell = a * np.diag(cell[:2])
if isinstance(cell, tuple):
cell = np.diag(cell)
slab.set_positions(positions.reshape((-1, 3)))
slab.set_cell([a * v * n for v, n in zip(cell, size)], scale_atoms=True)
if not periodic:
slab.cell[2] = 0.0
if vacuum is not None:
slab.center(vacuum, axis=2)
if 'adsorbate_info' not in slab.info:
slab.info.update({'adsorbate_info': {}})
slab.info['adsorbate_info']['cell'] = surface_cell
slab.info['adsorbate_info']['sites'] = sites
return slab
[docs]def fcc211(symbol, size, a=None, vacuum=None, orthogonal=True):
"""FCC(211) surface.
Does not currently support special adsorption sites.
Currently only implemented for *orthogonal=True* with size specified
as (i, j, k), where i, j, and k are number of atoms in each direction.
i must be divisible by 3 to accommodate the step width.
"""
if not orthogonal:
raise NotImplementedError('Only implemented for orthogonal '
'unit cells.')
if size[0] % 3 != 0:
raise NotImplementedError('First dimension of size must be '
'divisible by 3.')
atoms = FaceCenteredCubic(symbol,
directions=[[1, -1, -1],
[0, 2, -2],
[2, 1, 1]],
miller=(None, None, (2, 1, 1)),
latticeconstant=a,
size=(1, 1, 1),
pbc=True)
z = (size[2] + 1) // 2
atoms = atoms.repeat((size[0] // 3, size[1], z))
if size[2] % 2: # Odd: remove bottom layer and shrink cell.
remove_list = [atom.index for atom in atoms
if atom.z < atoms[1].z]
del atoms[remove_list]
dz = atoms[0].z
atoms.translate((0., 0., -dz))
atoms.cell[2][2] -= dz
atoms.cell[2] = 0.0
atoms.pbc[2] = False
if vacuum:
atoms.center(vacuum, axis=2)
# Renumber systematically from top down.
orders = [(atom.index, round(atom.x, 3), round(atom.y, 3),
-round(atom.z, 3), atom.index) for atom in atoms]
orders.sort(key=itemgetter(3, 1, 2))
newatoms = atoms.copy()
for index, order in enumerate(orders):
newatoms[index].position = atoms[order[0]].position.copy()
# Add empty 'sites' dictionary for consistency with other functions
newatoms.info['adsorbate_info'] = {'sites': {}}
return newatoms
[docs]def mx2(formula='MoS2', kind='2H', a=3.18, thickness=3.19,
size=(1, 1, 1), vacuum=None):
"""Create three-layer 2D materials with hexagonal structure.
For metal dichalcogenites, etc.
The kind argument accepts '2H', which gives a mirror plane symmetry
and '1T', which gives an inversion symmetry."""
if kind == '2H':
basis = [(0, 0, 0),
(2 / 3, 1 / 3, 0.5 * thickness),
(2 / 3, 1 / 3, -0.5 * thickness)]
elif kind == '1T':
basis = [(0, 0, 0),
(2 / 3, 1 / 3, 0.5 * thickness),
(1 / 3, 2 / 3, -0.5 * thickness)]
else:
raise ValueError('Structure not recognized:', kind)
cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
atoms.set_scaled_positions(basis)
if vacuum is not None:
atoms.center(vacuum, axis=2)
atoms = atoms.repeat(size)
return atoms
def graphene(formula='C2', a=2.460, size=(1, 1, 1), vacuum=None):
"""Create a graphene monolayer structure."""
cell = [[a, 0, 0], [-a / 2, a * 3**0.5 / 2, 0], [0, 0, 0]]
basis = [[0, 0, 0], [2 / 3, 1 / 3, 0]]
atoms = Atoms(formula, cell=cell, pbc=(1, 1, 0))
atoms.set_scaled_positions(basis)
if vacuum is not None:
atoms.center(vacuum, axis=2)
atoms = atoms.repeat(size)
return atoms
def _all_surface_functions():
# Convenient for debugging.
d = {}
for func in [fcc100, fcc110, bcc100, bcc110, bcc111, fcc111, hcp0001,
hcp10m10, diamond100, diamond111, fcc111, mx2, graphene]:
d[func.__name__] = func
return d