1from typing import List, Optional 

2 

3import numpy as np 

4 

5from ase.data import atomic_numbers as ref_atomic_numbers 

6from ase.spacegroup import Spacegroup 

7from ase.cluster.base import ClusterBase 

8from ase.cluster.cluster import Cluster 

9 

10 

11class ClusterFactory(ClusterBase): 

12 directions = [[1, 0, 0], 

13 [0, 1, 0], 

14 [0, 0, 1]] 

15 

16 atomic_basis = np.array([[0., 0., 0.]]) 

17 

18 element_basis: Optional[List[int]] = None 

19 

20 # Make it possible to change the class of the object returned. 

21 Cluster = Cluster 

22 

23 def __call__(self, symbols, surfaces, layers, latticeconstant=None, 

24 center=None, vacuum=0.0, debug=0): 

25 self.debug = debug 

26 

27 # Interpret symbol 

28 self.set_atomic_numbers(symbols) 

29 

30 # Interpret lattice constant 

31 if latticeconstant is None: 

32 if self.element_basis is None: 

33 self.lattice_constant = self.get_lattice_constant() 

34 else: 

35 raise ValueError( 

36 "A lattice constant must be specified for a compound") 

37 else: 

38 self.lattice_constant = latticeconstant 

39 

40 self.set_basis() 

41 

42 if self.debug: 

43 print("Lattice constant(s):", self.lattice_constant) 

44 print("Lattice basis:\n", self.lattice_basis) 

45 print("Resiprocal basis:\n", self.resiproc_basis) 

46 print("Atomic basis:\n", self.atomic_basis) 

47 

48 self.set_surfaces_layers(surfaces, layers) 

49 self.set_lattice_size(center) 

50 

51 if self.debug: 

52 print("Center position:", self.center.round(2)) 

53 print("Base lattice size:", self.size) 

54 

55 cluster = self.make_cluster(vacuum) 

56 cluster.symmetry = self.xtal_name 

57 cluster.surfaces = self.surfaces.copy() 

58 cluster.lattice_basis = self.lattice_basis.copy() 

59 cluster.atomic_basis = self.atomic_basis.copy() 

60 cluster.resiproc_basis = self.resiproc_basis.copy() 

61 return cluster 

62 

63 def make_cluster(self, vacuum): 

64 # Make the base crystal by repeating the unit cell 

65 size = np.array(self.size) 

66 translations = np.zeros((size.prod(), 3)) 

67 for h in range(size[0]): 

68 for k in range(size[1]): 

69 for l in range(size[2]): 

70 i = h * (size[1] * size[2]) + k * size[2] + l 

71 translations[i] = np.dot([h, k, l], self.lattice_basis) 

72 

73 atomic_basis = np.dot(self.atomic_basis, self.lattice_basis) 

74 positions = np.zeros((len(translations) * len(atomic_basis), 3)) 

75 numbers = np.zeros(len(positions)) 

76 n = len(atomic_basis) 

77 for i, trans in enumerate(translations): 

78 positions[n * i:n * (i + 1)] = atomic_basis + trans 

79 numbers[n * i:n * (i + 1)] = self.atomic_numbers 

80 

81 # Remove all atoms that is outside the defined surfaces 

82 for s, l in zip(self.surfaces, self.layers): 

83 n = self.miller_to_direction(s) 

84 rmax = self.get_layer_distance(s, l + 0.1) 

85 

86 r = np.dot(positions - self.center, n) 

87 mask = np.less(r, rmax) 

88 

89 if self.debug > 1: 

90 print("Cutting %s at %i layers ~ %.3f A" % (s, l, rmax)) 

91 

92 positions = positions[mask] 

93 numbers = numbers[mask] 

94 

95 atoms = self.Cluster(symbols=numbers, positions=positions) 

96 

97 atoms.cell = (1, 1, 1) # XXX ugly hack to center around zero 

98 atoms.center(about=(0, 0, 0)) 

99 atoms.cell[:] = 0 

100 return atoms 

101 

102 def set_atomic_numbers(self, symbols): 

103 "Extract atomic number from element" 

104 # The types that can be elements: integers and strings 

105 atomic_numbers = [] 

106 if self.element_basis is None: 

107 if isinstance(symbols, str): 

108 atomic_numbers.append(ref_atomic_numbers[symbols]) 

109 elif isinstance(symbols, int): 

110 atomic_numbers.append(symbols) 

111 else: 

112 raise TypeError("The symbol argument must be a " + 

113 "string or an atomic number.") 

114 element_basis = [0] * len(self.atomic_basis) 

115 else: 

116 if isinstance(symbols, (list, tuple)): 

117 nsymbols = len(symbols) 

118 else: 

119 nsymbols = 0 

120 

121 nelement_basis = max(self.element_basis) + 1 

122 if nsymbols != nelement_basis: 

123 raise TypeError("The symbol argument must be a sequence " + 

124 "of length %d" % (nelement_basis,) + 

125 " (one for each kind of lattice position") 

126 

127 for s in symbols: 

128 if isinstance(s, str): 

129 atomic_numbers.append(ref_atomic_numbers[s]) 

130 elif isinstance(s, int): 

131 atomic_numbers.append(s) 

132 else: 

133 raise TypeError("The symbol argument must be a " + 

134 "string or an atomic number.") 

135 element_basis = self.element_basis 

136 

137 self.atomic_numbers = [atomic_numbers[n] for n in element_basis] 

138 assert len(self.atomic_numbers) == len(self.atomic_basis) 

139 

140 def set_lattice_size(self, center): 

141 if center is None: 

142 offset = np.zeros(3) 

143 else: 

144 offset = np.array(center) 

145 if (offset > 1.0).any() or (offset < 0.0).any(): 

146 raise ValueError("Center offset must lie within the lattice unit \ 

147 cell.") 

148 

149 max = np.ones(3) 

150 min = -np.ones(3) 

151 v = np.linalg.inv(self.lattice_basis.T) 

152 for s, l in zip(self.surfaces, self.layers): 

153 n = self.miller_to_direction(s) * self.get_layer_distance(s, l) 

154 k = np.round(np.dot(v, n), 2) 

155 for i in range(3): 

156 if k[i] > 0.0: 

157 k[i] = np.ceil(k[i]) 

158 elif k[i] < 0.0: 

159 k[i] = np.floor(k[i]) 

160 

161 if self.debug > 1: 

162 print( 

163 "Spaning %i layers in %s in lattice basis ~ %s" % 

164 (l, s, k)) 

165 

166 max[k > max] = k[k > max] 

167 min[k < min] = k[k < min] 

168 

169 self.center = np.dot(offset - min, self.lattice_basis) 

170 self.size = (max - min + np.ones(3)).astype(int) 

171 

172 def set_surfaces_layers(self, surfaces, layers): 

173 if len(surfaces) != len(layers): 

174 raise ValueError( 

175 "Improper size of surface and layer arrays: %i != %i" 

176 % (len(surfaces), len(layers))) 

177 

178 sg = Spacegroup(self.spacegroup) 

179 surfaces = np.array(surfaces) 

180 layers = np.array(layers) 

181 

182 for i, s in enumerate(surfaces): 

183 s = reduce_miller(s) 

184 surfaces[i] = s 

185 

186 surfaces_full = surfaces.copy() 

187 layers_full = layers.copy() 

188 

189 for s, l in zip(surfaces, layers): 

190 equivalent_surfaces = sg.equivalent_reflections(s.reshape(-1, 3)) 

191 

192 for es in equivalent_surfaces: 

193 # If the equivalent surface (es) is not in the surface list, 

194 # then append it. 

195 if not np.equal(es, surfaces_full).all(axis=1).any(): 

196 surfaces_full = np.append( 

197 surfaces_full, es.reshape(1, 3), axis=0) 

198 layers_full = np.append(layers_full, l) 

199 

200 self.surfaces = surfaces_full.copy() 

201 self.layers = layers_full.copy() 

202 

203 def get_resiproc_basis(self, basis): 

204 """Returns the resiprocal basis to a given lattice (crystal) basis""" 

205 k = 1 / np.dot(basis[0], cross(basis[1], basis[2])) 

206 

207 # The same as the inversed basis matrix transposed 

208 return k * np.array([cross(basis[1], basis[2]), 

209 cross(basis[2], basis[0]), 

210 cross(basis[0], basis[1])]) 

211 

212# Helping functions 

213 

214 

215def cross(a, b): 

216 """The cross product of two vectors.""" 

217 return np.array([a[1] * b[2] - b[1] * a[2], 

218 a[2] * b[0] - b[2] * a[0], 

219 a[0] * b[1] - b[0] * a[1]]) 

220 

221 

222def GCD(a, b): 

223 """Greatest Common Divisor of a and b.""" 

224 # print "--" 

225 while a != 0: 

226 # print a,b,">", 

227 a, b = b % a, a 

228 # print a,b 

229 return b 

230 

231 

232def reduce_miller(hkl): 

233 """Reduce Miller index to the lowest equivalent integers.""" 

234 hkl = np.array(hkl) 

235 old = hkl.copy() 

236 

237 d = GCD(GCD(hkl[0], hkl[1]), hkl[2]) 

238 while d != 1: 

239 hkl = hkl // d 

240 d = GCD(GCD(hkl[0], hkl[1]), hkl[2]) 

241 

242 if np.dot(old, hkl) > 0: 

243 return hkl 

244 else: 

245 return -hkl