Coverage for /builds/debichem-team/python-ase/ase/build/general_surface.py: 83.61%

1from math import gcd 

2 

3import numpy as np 

4from numpy.linalg import norm, solve 

5 

6from ase.build import bulk 

7from ase.build.surface import create_tags 

8 

9 

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

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

12 

13 lattice: Atoms object or str 

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

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

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

17 correct bulk lattice will be generated automatically. 

18 indices: sequence of three int 

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

20 layers: int 

21 Number of equivalent layers of the slab. 

22 vacuum: float 

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

24 periodic: bool 

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

26 """ 

27 

28 indices = np.asarray(indices) 

29 

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

31 raise ValueError(f'{indices} is an invalid surface type') 

32 

33 if isinstance(lattice, str): 

34 lattice = bulk(lattice, cubic=True) 

35 

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

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

38 

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

40 if not h0: 

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

42 if not k0: 

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

44 if not l0: 

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

46 else: 

47 p, q = ext_gcd(k, l) 

48 a1, a2, a3 = lattice.cell 

49 

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

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

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

53 l * a2 - k * a3) 

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

55 l * a2 - k * a3) 

56 

57 if abs(k2) > tol: 

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

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

60 

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

62 

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

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

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

66 

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

68 if vacuum is not None: 

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

70 return surf 

71 

72 

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

74 surf = lattice.copy() 

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

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

77 surf.set_scaled_positions(scaled) 

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

79 surf *= (1, 1, layers) 

80 surf.set_tags(create_tags((1, len(lattice), layers))) 

81 

82 a1, a2, a3 = surf.cell 

83 surf.set_cell([a1, a2, 

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

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

86 

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

88 # and z perpendicular to the surface: 

89 a1, a2, a3 = surf.cell 

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

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

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

93 (0, 0, norm(a3))], 

94 scale_atoms=True) 

95 

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

97 

98 # Move atoms into the unit cell: 

99 scaled = surf.get_scaled_positions() 

100 scaled[:, :2] %= 1 

101 surf.set_scaled_positions(scaled) 

102 

103 if not periodic: 

104 surf.cell[2] = 0.0 

105 

106 return surf 

107 

108 

109def ext_gcd(a, b): 

110 if b == 0: 

111 return 1, 0 

112 elif a % b == 0: 

113 return 0, 1 

114 else: 

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

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