1"""Helper functions for creating supercells.""" 

2 

3import numpy as np 

4 

5from ase import Atoms 

6 

7 

8class SupercellError(Exception): 

9 """Use if construction of supercell fails""" 

10 

11 

12def get_deviation_from_optimal_cell_shape(cell, target_shape="sc", norm=None): 

13 r""" 

14 Calculates the deviation of the given cell metric from the ideal 

15 cell metric defining a certain shape. Specifically, the function 

16 evaluates the expression `\Delta = || Q \mathbf{h} - 

17 \mathbf{h}_{target}||_2`, where `\mathbf{h}` is the input 

18 metric (*cell*) and `Q` is a normalization factor (*norm*) 

19 while the target metric `\mathbf{h}_{target}` (via 

20 *target_shape*) represent simple cubic ('sc') or face-centered 

21 cubic ('fcc') cell shapes. 

22 

23 Parameters: 

24 

25 cell: 2D array of floats 

26 Metric given as a (3x3 matrix) of the input structure. 

27 target_shape: str 

28 Desired supercell shape. Can be 'sc' for simple cubic or 

29 'fcc' for face-centered cubic. 

30 norm: float 

31 Specify the normalization factor. This is useful to avoid 

32 recomputing the normalization factor when computing the 

33 deviation for a series of P matrices. 

34 

35 """ 

36 

37 if target_shape in ["sc", "simple-cubic"]: 

38 target_metric = np.eye(3) 

39 elif target_shape in ["fcc", "face-centered cubic"]: 

40 target_metric = 0.5 * np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) 

41 if not norm: 

42 norm = (np.linalg.det(cell) / np.linalg.det(target_metric)) ** ( 

43 -1.0 / 3 

44 ) 

45 return np.linalg.norm(norm * cell - target_metric) 

46 

47 

48def find_optimal_cell_shape( 

49 cell, 

50 target_size, 

51 target_shape, 

52 lower_limit=-2, 

53 upper_limit=2, 

54 verbose=False, 

55): 

56 """Returns the transformation matrix that produces a supercell 

57 corresponding to *target_size* unit cells with metric *cell* that 

58 most closely approximates the shape defined by *target_shape*. 

59 

60 Parameters: 

61 

62 cell: 2D array of floats 

63 Metric given as a (3x3 matrix) of the input structure. 

64 target_size: integer 

65 Size of desired super cell in number of unit cells. 

66 target_shape: str 

67 Desired supercell shape. Can be 'sc' for simple cubic or 

68 'fcc' for face-centered cubic. 

69 lower_limit: int 

70 Lower limit of search range. 

71 upper_limit: int 

72 Upper limit of search range. 

73 verbose: bool 

74 Set to True to obtain additional information regarding 

75 construction of transformation matrix. 

76 

77 """ 

78 

79 # Set up target metric 

80 if target_shape in ["sc", "simple-cubic"]: 

81 target_metric = np.eye(3) 

82 elif target_shape in ["fcc", "face-centered cubic"]: 

83 target_metric = 0.5 * np.array( 

84 [[0, 1, 1], [1, 0, 1], [1, 1, 0]], dtype=float 

85 ) 

86 if verbose: 

87 print("target metric (h_target):") 

88 print(target_metric) 

89 

90 # Normalize cell metric to reduce computation time during looping 

91 norm = ( 

92 target_size * np.linalg.det(cell) / np.linalg.det(target_metric) 

93 ) ** (-1.0 / 3) 

94 norm_cell = norm * cell 

95 if verbose: 

96 print("normalization factor (Q): %g" % norm) 

97 

98 # Approximate initial P matrix 

99 ideal_P = np.dot(target_metric, np.linalg.inv(norm_cell)) 

100 if verbose: 

101 print("idealized transformation matrix:") 

102 print(ideal_P) 

103 starting_P = np.array(np.around(ideal_P, 0), dtype=int) 

104 if verbose: 

105 print("closest integer transformation matrix (P_0):") 

106 print(starting_P) 

107 

108 # Prepare run. 

109 from itertools import product 

110 

111 best_score = 1e6 

112 optimal_P = None 

113 for dP in product(range(lower_limit, upper_limit + 1), repeat=9): 

114 dP = np.array(dP, dtype=int).reshape(3, 3) 

115 P = starting_P + dP 

116 if int(np.around(np.linalg.det(P), 0)) != target_size: 

117 continue 

118 score = get_deviation_from_optimal_cell_shape( 

119 np.dot(P, norm_cell), target_shape=target_shape, norm=1.0 

120 ) 

121 if score < best_score: 

122 best_score = score 

123 optimal_P = P 

124 

125 if optimal_P is None: 

126 print("Failed to find a transformation matrix.") 

127 return None 

128 

129 # Finalize. 

130 if verbose: 

131 print("smallest score (|Q P h_p - h_target|_2): %f" % best_score) 

132 print("optimal transformation matrix (P_opt):") 

133 print(optimal_P) 

134 print("supercell metric:") 

135 print(np.round(np.dot(optimal_P, cell), 4)) 

136 print( 

137 "determinant of optimal transformation matrix: %g" 

138 % np.linalg.det(optimal_P) 

139 ) 

140 return optimal_P 

141 

142 

143def make_supercell(prim, P, wrap=True, tol=1e-5): 

144 r"""Generate a supercell by applying a general transformation (*P*) to 

145 the input configuration (*prim*). 

146 

147 The transformation is described by a 3x3 integer matrix 

148 `\mathbf{P}`. Specifically, the new cell metric 

149 `\mathbf{h}` is given in terms of the metric of the input 

150 configuration `\mathbf{h}_p` by `\mathbf{P h}_p = 

151 \mathbf{h}`. 

152 

153 Parameters: 

154 

155 prim: ASE Atoms object 

156 Input configuration. 

157 P: 3x3 integer matrix 

158 Transformation matrix `\mathbf{P}`. 

159 wrap: bool 

160 wrap in the end 

161 tol: float 

162 tolerance for wrapping 

163 """ 

164 

165 supercell_matrix = P 

166 supercell = clean_matrix(supercell_matrix @ prim.cell) 

167 

168 # cartesian lattice points 

169 lattice_points_frac = lattice_points_in_supercell(supercell_matrix) 

170 lattice_points = np.dot(lattice_points_frac, supercell) 

171 

172 superatoms = Atoms(cell=supercell, pbc=prim.pbc) 

173 

174 for lp in lattice_points: 

175 shifted_atoms = prim.copy() 

176 shifted_atoms.positions += lp 

177 superatoms.extend(shifted_atoms) 

178 

179 # check number of atoms is correct 

180 n_target = int(np.round(np.linalg.det(supercell_matrix) * len(prim))) 

181 if n_target != len(superatoms): 

182 msg = "Number of atoms in supercell: {}, expected: {}".format( 

183 n_target, len(superatoms) 

184 ) 

185 raise SupercellError(msg) 

186 

187 if wrap: 

188 superatoms.wrap(eps=tol) 

189 

190 return superatoms 

191 

192 

193def lattice_points_in_supercell(supercell_matrix): 

194 """Find all lattice points contained in a supercell. 

195 

196 Adapted from pymatgen, which is available under MIT license: 

197 The MIT License (MIT) Copyright (c) 2011-2012 MIT & The Regents of the 

198 University of California, through Lawrence Berkeley National Laboratory 

199 """ 

200 

201 diagonals = np.array( 

202 [ 

203 [0, 0, 0], 

204 [0, 0, 1], 

205 [0, 1, 0], 

206 [0, 1, 1], 

207 [1, 0, 0], 

208 [1, 0, 1], 

209 [1, 1, 0], 

210 [1, 1, 1], 

211 ] 

212 ) 

213 d_points = np.dot(diagonals, supercell_matrix) 

214 

215 mins = np.min(d_points, axis=0) 

216 maxes = np.max(d_points, axis=0) + 1 

217 

218 ar = np.arange(mins[0], maxes[0])[:, None] * np.array([1, 0, 0])[None, :] 

219 br = np.arange(mins[1], maxes[1])[:, None] * np.array([0, 1, 0])[None, :] 

220 cr = np.arange(mins[2], maxes[2])[:, None] * np.array([0, 0, 1])[None, :] 

221 

222 all_points = ar[:, None, None] + br[None, :, None] + cr[None, None, :] 

223 all_points = all_points.reshape((-1, 3)) 

224 

225 frac_points = np.dot(all_points, np.linalg.inv(supercell_matrix)) 

226 

227 tvects = frac_points[ 

228 np.all(frac_points < 1 - 1e-10, axis=1) 

229 & np.all(frac_points >= -1e-10, axis=1) 

230 ] 

231 assert len(tvects) == round(abs(np.linalg.det(supercell_matrix))) 

232 return tvects 

233 

234 

235def clean_matrix(matrix, eps=1e-12): 

236 """ clean from small values""" 

237 matrix = np.array(matrix) 

238 for ij in np.ndindex(matrix.shape): 

239 if abs(matrix[ij]) < eps: 

240 matrix[ij] = 0 

241 return matrix