r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1from math import gcd

2import numpy as np

3from numpy.linalg import norm, solve

5from ase.build import bulk

8def surface(lattice, indices, layers, vacuum=None, tol=1e-10, periodic=False):

9 """Create surface from a given lattice and Miller indices.

11 lattice: Atoms object or str

12 Bulk lattice structure of alloy or pure metal. Note that the

13 unit-cell must be the conventional cell - not the primitive cell.

14 One can also give the chemical symbol as a string, in which case the

15 correct bulk lattice will be generated automatically.

16 indices: sequence of three int

17 Surface normal in Miller indices (h,k,l).

18 layers: int

19 Number of equivalent layers of the slab.

20 vacuum: float

21 Amount of vacuum added on both sides of the slab.

22 periodic: bool

23 Whether the surface is periodic in the normal to the surface

24 """

26 indices = np.asarray(indices)

28 if indices.shape != (3,) or not indices.any() or indices.dtype != int:

29 raise ValueError('%s is an invalid surface type' % indices)

31 if isinstance(lattice, str):

32 lattice = bulk(lattice, cubic=True)

34 h, k, l = indices # noqa (E741, the variable l)

35 h0, k0, l0 = (indices == 0)

37 if h0 and k0 or h0 and l0 or k0 and l0: # if two indices are zero

38 if not h0:

39 c1, c2, c3 = [(0, 1, 0), (0, 0, 1), (1, 0, 0)]

40 if not k0:

41 c1, c2, c3 = [(0, 0, 1), (1, 0, 0), (0, 1, 0)]

42 if not l0:

43 c1, c2, c3 = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]

44 else:

45 p, q = ext_gcd(k, l)

46 a1, a2, a3 = lattice.cell

48 # constants describing the dot product of basis c1 and c2:

49 # dot(c1,c2) = k1+i*k2, i in Z

50 k1 = np.dot(p * (k * a1 - h * a2) + q * (l * a1 - h * a3),

51 l * a2 - k * a3)

52 k2 = np.dot(l * (k * a1 - h * a2) - k * (l * a1 - h * a3),

53 l * a2 - k * a3)

55 if abs(k2) > tol:

56 i = -int(round(k1 / k2)) # i corresponding to the optimal basis

57 p, q = p + i * l, q - i * k

59 a, b = ext_gcd(p * k + q * l, h)

61 c1 = (p * k + q * l, -p * h, -q * h)

62 c2 = np.array((0, l, -k)) // abs(gcd(l, k))

63 c3 = (b, a * p, a * q)

65 surf = build(lattice, np.array([c1, c2, c3]), layers, tol, periodic)

66 if vacuum is not None:

67 surf.center(vacuum=vacuum, axis=2)

68 return surf

71def build(lattice, basis, layers, tol, periodic):

72 surf = lattice.copy()

73 scaled = solve(basis.T, surf.get_scaled_positions().T).T

74 scaled -= np.floor(scaled + tol)

75 surf.set_scaled_positions(scaled)

76 surf.set_cell(np.dot(basis, surf.cell), scale_atoms=True)

77 surf *= (1, 1, layers)

79 a1, a2, a3 = surf.cell

80 surf.set_cell([a1, a2,

81 np.cross(a1, a2) * np.dot(a3, np.cross(a1, a2)) /

82 norm(np.cross(a1, a2))**2])

84 # Change unit cell to have the x-axis parallel with a surface vector

85 # and z perpendicular to the surface:

86 a1, a2, a3 = surf.cell

87 surf.set_cell([(norm(a1), 0, 0),

88 (np.dot(a1, a2) / norm(a1),

89 np.sqrt(norm(a2)**2 - (np.dot(a1, a2) / norm(a1))**2), 0),

90 (0, 0, norm(a3))],

91 scale_atoms=True)

93 surf.pbc = (True, True, periodic)

95 # Move atoms into the unit cell:

96 scaled = surf.get_scaled_positions()

97 scaled[:, :2] %= 1

98 surf.set_scaled_positions(scaled)

100 if not periodic:

101 surf.cell[2] = 0.0

103 return surf

106def ext_gcd(a, b):

107 if b == 0:

108 return 1, 0

109 elif a % b == 0:

110 return 0, 1

111 else:

112 x, y = ext_gcd(b, a % b)

113 return y, x - y * (a // b)