1from math import log10, atan2, cos, sin 

2from ase.build import hcp0001, fcc111, bcc111 

3import numpy as np 

4 

5 

6def hcp0001_root(symbol, root, size, a=None, c=None, 

7 vacuum=None, orthogonal=False): 

8 """HCP(0001) surface maniupulated to have a x unit side length 

9 of *root* before repeating. This also results in *root* number 

10 of repetitions of the cell. 

11 

12 

13 The first 20 valid roots for nonorthogonal are... 

14 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25, 

15 27, 28, 31, 36, 37, 39, 43, 48, 49""" 

16 atoms = hcp0001(symbol=symbol, size=(1, 1, size[2]), 

17 a=a, c=c, vacuum=vacuum, orthogonal=orthogonal) 

18 atoms = root_surface(atoms, root) 

19 atoms *= (size[0], size[1], 1) 

20 return atoms 

21 

22 

23def fcc111_root(symbol, root, size, a=None, 

24 vacuum=None, orthogonal=False): 

25 """FCC(111) surface maniupulated to have a x unit side length 

26 of *root* before repeating. This also results in *root* number 

27 of repetitions of the cell. 

28 

29 The first 20 valid roots for nonorthogonal are... 

30 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25, 27, 

31 28, 31, 36, 37, 39, 43, 48, 49""" 

32 atoms = fcc111(symbol=symbol, size=(1, 1, size[2]), 

33 a=a, vacuum=vacuum, orthogonal=orthogonal) 

34 atoms = root_surface(atoms, root) 

35 atoms *= (size[0], size[1], 1) 

36 return atoms 

37 

38 

39def bcc111_root(symbol, root, size, a=None, 

40 vacuum=None, orthogonal=False): 

41 """BCC(111) surface maniupulated to have a x unit side length 

42 of *root* before repeating. This also results in *root* number 

43 of repetitions of the cell. 

44 

45 

46 The first 20 valid roots for nonorthogonal are... 

47 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25, 

48 27, 28, 31, 36, 37, 39, 43, 48, 49""" 

49 atoms = bcc111(symbol=symbol, size=(1, 1, size[2]), 

50 a=a, vacuum=vacuum, orthogonal=orthogonal) 

51 atoms = root_surface(atoms, root) 

52 atoms *= (size[0], size[1], 1) 

53 return atoms 

54 

55 

56def point_in_cell_2d(point, cell, eps=1e-8): 

57 """This function takes a 2D slice of the cell in the XY plane and calculates 

58 if a point should lie in it. This is used as a more accurate method of 

59 ensuring we find all of the correct cell repetitions in the root surface 

60 code. The Z axis is totally ignored but for most uses this should be fine. 

61 """ 

62 # Define area of a triangle 

63 def tri_area(t1, t2, t3): 

64 t1x, t1y = t1[0:2] 

65 t2x, t2y = t2[0:2] 

66 t3x, t3y = t3[0:2] 

67 return abs(t1x * (t2y - t3y) + t2x * (t3y - t1y) + t3x * (t1y - t2y)) / 2 

68 

69 # c0, c1, c2, c3 define a parallelogram 

70 c0 = (0, 0) 

71 c1 = cell[0, 0:2] 

72 c2 = cell[1, 0:2] 

73 c3 = c1 + c2 

74 

75 # Get area of parallelogram 

76 cA = tri_area(c0, c1, c2) + tri_area(c1, c2, c3) 

77 

78 # Get area of triangles formed from adjacent vertices of parallelogram and 

79 # point in question. 

80 pA = tri_area(point, c0, c1) + tri_area(point, c1, c2) + tri_area(point, c2, c3) + tri_area(point, c3, c0) 

81 

82 # If combined area of triangles from point is larger than area of 

83 # parallelogram, point is not inside parallelogram. 

84 return pA <= cA + eps 

85 

86 

87def _root_cell_normalization(primitive_slab): 

88 """Returns the scaling factor for x axis and cell normalized by that factor""" 

89 

90 xscale = np.linalg.norm(primitive_slab.cell[0, 0:2]) 

91 cell_vectors = primitive_slab.cell[0:2, 0:2] / xscale 

92 return xscale, cell_vectors 

93 

94 

95def _root_surface_analysis(primitive_slab, root, eps=1e-8): 

96 """A tool to analyze a slab and look for valid roots that exist, up to 

97 the given root. This is useful for generating all possible cells 

98 without prior knowledge. 

99 

100 *primitive slab* is the primitive cell to analyze. 

101 

102 *root* is the desired root to find, and all below. 

103 

104 This is the internal function which gives extra data to root_surface. 

105 """ 

106 

107 # Setup parameters for cell searching 

108 logeps = int(-log10(eps)) 

109 xscale, cell_vectors = _root_cell_normalization(primitive_slab) 

110 

111 # Allocate grid for cell search search 

112 points = np.indices((root + 1, root + 1)).T.reshape(-1, 2) 

113 

114 # Find points corresponding to full cells 

115 cell_points = [cell_vectors[0] * x + cell_vectors[1] * y for x, y in points] 

116 

117 # Find point close to the desired cell (floating point error possible) 

118 roots = np.around(np.linalg.norm(cell_points, axis=1)**2, logeps) 

119 

120 valid_roots = np.nonzero(roots == root)[0] 

121 if len(valid_roots) == 0: 

122 raise ValueError("Invalid root {} for cell {}".format(root, cell_vectors)) 

123 int_roots = np.array([int(this_root) for this_root in roots 

124 if this_root.is_integer() and this_root <= root]) 

125 return cell_points, cell_points[np.nonzero(roots == root)[0][0]], set(int_roots[1:]) 

126 

127 

128def root_surface_analysis(primitive_slab, root, eps=1e-8): 

129 """A tool to analyze a slab and look for valid roots that exist, up to 

130 the given root. This is useful for generating all possible cells 

131 without prior knowledge. 

132 

133 *primitive slab* is the primitive cell to analyze. 

134 

135 *root* is the desired root to find, and all below.""" 

136 return _root_surface_analysis(primitive_slab=primitive_slab, root=root, eps=eps)[2] 

137 

138 

139def root_surface(primitive_slab, root, eps=1e-8): 

140 """Creates a cell from a primitive cell that repeats along the x and y 

141 axis in a way consisent with the primitive cell, that has been cut 

142 to have a side length of *root*. 

143 

144 *primitive cell* should be a primitive 2d cell of your slab, repeated 

145 as needed in the z direction. 

146 

147 *root* should be determined using an analysis tool such as the 

148 root_surface_analysis function, or prior knowledge. It should always 

149 be a whole number as it represents the number of repetitions.""" 

150 

151 atoms = primitive_slab.copy() 

152 

153 xscale, cell_vectors = _root_cell_normalization(primitive_slab) 

154 

155 # Do root surface analysis 

156 cell_points, root_point, roots = _root_surface_analysis(primitive_slab, root, eps=eps) 

157 

158 # Find new cell 

159 root_angle = -atan2(root_point[1], root_point[0]) 

160 root_rotation = [[cos(root_angle), -sin(root_angle)], 

161 [sin(root_angle), cos(root_angle)]] 

162 root_scale = np.linalg.norm(root_point) 

163 

164 cell = np.array([np.dot(x, root_rotation) * root_scale for x in cell_vectors]) 

165 

166 # Find all cell centers within the cell 

167 shift = cell_vectors.sum(axis=0) / 2 

168 cell_points = [point for point in cell_points if point_in_cell_2d(point+shift, cell, eps=eps)] 

169 

170 # Setup new cell 

171 atoms.rotate(root_angle, v="z") 

172 atoms *= (root, root, 1) 

173 atoms.cell[0:2, 0:2] = cell * xscale 

174 atoms.center() 

175 

176 # Remove all extra atoms 

177 del atoms[[atom.index for atom in atoms if not point_in_cell_2d(atom.position, atoms.cell, eps=eps)]] 

178 

179 # Rotate cell back to original orientation 

180 standard_rotation = [[cos(-root_angle), -sin(-root_angle), 0], 

181 [sin(-root_angle), cos(-root_angle), 0], 

182 [0, 0, 1]] 

183 

184 new_cell = np.array([np.dot(x, standard_rotation) for x in atoms.cell]) 

185 new_positions = np.array([np.dot(x, standard_rotation) 

186 for x in atoms.positions]) 

187 

188 atoms.cell = new_cell 

189 atoms.positions = new_positions 

190 return atoms