Coverage for ase / cluster / factory.py: 80.00%

1# fmt: off 

2 

3 

4import numpy as np 

5 

6from ase.cluster.base import ClusterBase 

7from ase.cluster.cluster import Cluster 

8from ase.data import atomic_numbers as ref_atomic_numbers 

9from ase.spacegroup import Spacegroup 

10 

11 

12class ClusterFactory(ClusterBase): 

13 directions = [[1, 0, 0], 

14 [0, 1, 0], 

15 [0, 0, 1]] 

16 

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

18 

19 element_basis: list[int] | None = None 

20 

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

22 Cluster = Cluster 

23 

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

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

26 self.debug = debug 

27 

28 # Interpret symbol 

29 self.set_atomic_numbers(symbols) 

30 

31 # Interpret lattice constant 

32 if latticeconstant is None: 

33 if self.element_basis is None: 

34 self.lattice_constant = self.get_lattice_constant() 

35 else: 

36 raise ValueError( 

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

38 else: 

39 self.lattice_constant = latticeconstant 

40 

41 self.set_basis() 

42 

43 if self.debug: 

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

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

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

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

48 

49 self.set_surfaces_layers(surfaces, layers) 

50 self.set_lattice_size(center) 

51 

52 if self.debug: 

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

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

55 

56 cluster = self.make_cluster(vacuum) 

57 cluster.symmetry = self.xtal_name 

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

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

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

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

62 return cluster 

63 

64 def make_cluster(self, vacuum): 

65 # Make the base crystal by repeating the unit cell 

66 size = np.array(self.size) 

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

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

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

70 for l_ in range(size[2]): 

71 i = h * (size[1] * size[2]) + k * size[2] + l_ 

72 translations[i] = np.dot([h, k, l_], self.lattice_basis) 

73 

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

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

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

77 n = len(atomic_basis) 

78 for i, trans in enumerate(translations): 

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

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

81 

82 # Remove all atoms that is outside the defined surfaces 

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

84 n = self.miller_to_direction(s) 

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

86 

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

88 mask = np.less(r, rmax) 

89 

90 if self.debug > 1: 

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

92 

93 positions = positions[mask] 

94 numbers = numbers[mask] 

95 

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

97 

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

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

100 atoms.cell[:] = 0 

101 return atoms 

102 

103 def set_atomic_numbers(self, symbols): 

104 "Extract atomic number from element" 

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

106 atomic_numbers = [] 

107 if self.element_basis is None: 

108 if isinstance(symbols, str): 

109 atomic_numbers.append(ref_atomic_numbers[symbols]) 

110 elif isinstance(symbols, int): 

111 atomic_numbers.append(symbols) 

112 else: 

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

114 "string or an atomic number.") 

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

116 else: 

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

118 nsymbols = len(symbols) 

119 else: 

120 nsymbols = 0 

121 

122 nelement_basis = max(self.element_basis) + 1 

123 if nsymbols != nelement_basis: 

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

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

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

127 

128 for s in symbols: 

129 if isinstance(s, str): 

130 atomic_numbers.append(ref_atomic_numbers[s]) 

131 elif isinstance(s, int): 

132 atomic_numbers.append(s) 

133 else: 

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

135 "string or an atomic number.") 

136 element_basis = self.element_basis 

137 

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

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

140 

141 def set_lattice_size(self, center): 

142 if center is None: 

143 offset = np.zeros(3) 

144 else: 

145 offset = np.array(center) 

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

147 raise ValueError( 

148 "Center offset must lie within the lattice unit cell.") 

149 

150 max = np.ones(3) 

151 min = -np.ones(3) 

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

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

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

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

156 for i in range(3): 

157 if k[i] > 0.0: 

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

159 elif k[i] < 0.0: 

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

161 

162 if self.debug > 1: 

163 print( 

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

165 (l, s, k)) 

166 

167 max[k > max] = k[k > max] 

168 min[k < min] = k[k < min] 

169 

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

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

172 

173 def set_surfaces_layers(self, surfaces, layers): 

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

175 raise ValueError( 

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

177 % (len(surfaces), len(layers))) 

178 

179 sg = Spacegroup(self.spacegroup) 

180 surfaces = np.array(surfaces) 

181 layers = np.array(layers) 

182 

183 for i, s in enumerate(surfaces): 

184 s = reduce_miller(s) 

185 surfaces[i] = s 

186 

187 surfaces_full = surfaces.copy() 

188 layers_full = layers.copy() 

189 

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

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

192 

193 for es in equivalent_surfaces: 

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

195 # then append it. 

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

197 surfaces_full = np.append( 

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

199 layers_full = np.append(layers_full, l) 

200 

201 self.surfaces = surfaces_full.copy() 

202 self.layers = layers_full.copy() 

203 

204 def get_resiproc_basis(self, basis): 

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

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

207 

208 # The same as the inversed basis matrix transposed 

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

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

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

212 

213# Helping functions 

214 

215 

216def cross(a, b): 

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

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

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

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

221 

222 

223def GCD(a, b): 

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

225 # print "--" 

226 while a != 0: 

227 # print a,b,">", 

228 a, b = b % a, a 

229 # print a,b 

230 return b 

231 

232 

233def reduce_miller(hkl): 

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

235 hkl = np.array(hkl) 

236 old = hkl.copy() 

237 

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

239 while d != 1: 

240 hkl = hkl // d 

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

242 

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

244 return hkl 

245 else: 

246 return -hkl