Coverage for /builds/ase/ase/ase/cluster/factory.py: 80.13%

1# fmt: off 

2 

3from typing import List, Optional 

4 

5import numpy as np 

6 

7from ase.cluster.base import ClusterBase 

8from ase.cluster.cluster import Cluster 

9from ase.data import atomic_numbers as ref_atomic_numbers 

10from ase.spacegroup import Spacegroup 

11 

12 

13class ClusterFactory(ClusterBase): 

14 directions = [[1, 0, 0], 

15 [0, 1, 0], 

16 [0, 0, 1]] 

17 

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

19 

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

21 

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

23 Cluster = Cluster 

24 

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

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

27 self.debug = debug 

28 

29 # Interpret symbol 

30 self.set_atomic_numbers(symbols) 

31 

32 # Interpret lattice constant 

33 if latticeconstant is None: 

34 if self.element_basis is None: 

35 self.lattice_constant = self.get_lattice_constant() 

36 else: 

37 raise ValueError( 

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

39 else: 

40 self.lattice_constant = latticeconstant 

41 

42 self.set_basis() 

43 

44 if self.debug: 

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

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

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

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

49 

50 self.set_surfaces_layers(surfaces, layers) 

51 self.set_lattice_size(center) 

52 

53 if self.debug: 

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

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

56 

57 cluster = self.make_cluster(vacuum) 

58 cluster.symmetry = self.xtal_name 

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

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

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

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

63 return cluster 

64 

65 def make_cluster(self, vacuum): 

66 # Make the base crystal by repeating the unit cell 

67 size = np.array(self.size) 

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

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

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

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

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

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

74 

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

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

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

78 n = len(atomic_basis) 

79 for i, trans in enumerate(translations): 

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

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

82 

83 # Remove all atoms that is outside the defined surfaces 

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

85 n = self.miller_to_direction(s) 

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

87 

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

89 mask = np.less(r, rmax) 

90 

91 if self.debug > 1: 

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

93 

94 positions = positions[mask] 

95 numbers = numbers[mask] 

96 

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

98 

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

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

101 atoms.cell[:] = 0 

102 return atoms 

103 

104 def set_atomic_numbers(self, symbols): 

105 "Extract atomic number from element" 

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

107 atomic_numbers = [] 

108 if self.element_basis is None: 

109 if isinstance(symbols, str): 

110 atomic_numbers.append(ref_atomic_numbers[symbols]) 

111 elif isinstance(symbols, int): 

112 atomic_numbers.append(symbols) 

113 else: 

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

115 "string or an atomic number.") 

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

117 else: 

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

119 nsymbols = len(symbols) 

120 else: 

121 nsymbols = 0 

122 

123 nelement_basis = max(self.element_basis) + 1 

124 if nsymbols != nelement_basis: 

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

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

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

128 

129 for s in symbols: 

130 if isinstance(s, str): 

131 atomic_numbers.append(ref_atomic_numbers[s]) 

132 elif isinstance(s, int): 

133 atomic_numbers.append(s) 

134 else: 

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

136 "string or an atomic number.") 

137 element_basis = self.element_basis 

138 

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

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

141 

142 def set_lattice_size(self, center): 

143 if center is None: 

144 offset = np.zeros(3) 

145 else: 

146 offset = np.array(center) 

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

148 raise ValueError( 

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

150 

151 max = np.ones(3) 

152 min = -np.ones(3) 

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

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

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

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

157 for i in range(3): 

158 if k[i] > 0.0: 

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

160 elif k[i] < 0.0: 

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

162 

163 if self.debug > 1: 

164 print( 

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

166 (l, s, k)) 

167 

168 max[k > max] = k[k > max] 

169 min[k < min] = k[k < min] 

170 

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

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

173 

174 def set_surfaces_layers(self, surfaces, layers): 

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

176 raise ValueError( 

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

178 % (len(surfaces), len(layers))) 

179 

180 sg = Spacegroup(self.spacegroup) 

181 surfaces = np.array(surfaces) 

182 layers = np.array(layers) 

183 

184 for i, s in enumerate(surfaces): 

185 s = reduce_miller(s) 

186 surfaces[i] = s 

187 

188 surfaces_full = surfaces.copy() 

189 layers_full = layers.copy() 

190 

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

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

193 

194 for es in equivalent_surfaces: 

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

196 # then append it. 

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

198 surfaces_full = np.append( 

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

200 layers_full = np.append(layers_full, l) 

201 

202 self.surfaces = surfaces_full.copy() 

203 self.layers = layers_full.copy() 

204 

205 def get_resiproc_basis(self, basis): 

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

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

208 

209 # The same as the inversed basis matrix transposed 

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

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

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

213 

214# Helping functions 

215 

216 

217def cross(a, b): 

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

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

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

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

222 

223 

224def GCD(a, b): 

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

226 # print "--" 

227 while a != 0: 

228 # print a,b,">", 

229 a, b = b % a, a 

230 # print a,b 

231 return b 

232 

233 

234def reduce_miller(hkl): 

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

236 hkl = np.array(hkl) 

237 old = hkl.copy() 

238 

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

240 while d != 1: 

241 hkl = hkl // d 

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

243 

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

245 return hkl 

246 else: 

247 return -hkl