Coverage for hiphive/core/orbits.py: 97%

1""" 

2Contains the Orbit class which hold onformation about equivalent clusters. 

3""" 

4import numpy as np 

5import tarfile 

6import pickle 

7 

8from .utilities import BiMap 

9from .atoms import Atom 

10from ..input_output.logging_tools import logger 

11from ..input_output.read_write_files import add_items_to_tarfile_pickle, read_items_pickle 

12 

13 

14logger = logger.getChild('orbits') 

15 

16 

17# This is the interface accessible for cluster_space 

18def get_orbits(cluster_list, atom_list, rotation_matrices, translation_vectors, 

19 permutations, prim, symprec): 

20 """Generate a list of the orbits for the clusters in a supercell 

21 configuration. 

22 

23 This method requires as input a list of the clusters in a supercell 

24 configuration as well as a set of symmetry operations (rotations and 

25 translations). From this information it will generate a list of the orbits, 

26 i.e. the set of symmetry inequivalent clusters each associated with its 

27 respective set of equivalent clusters. 

28 

29 Parameters 

30 ---------- 

31 cluster_list : BiMap object 

32 a list of clusters 

33 atom_list : BiMap object 

34 a list of atoms in a supercell 

35 rotation_matrices : list of NumPy (3,3) arrays 

36 rotational symmetries to be imposed (e.g., from spglib) 

37 translation_vectors : list of NumPy (3) arrays 

38 translational symmetries to be imposed (e.g., from spglib) 

39 permutations : list of permutations 

40 lookup table for permutations 

41 prim : hiPhive Atoms object 

42 primitive structure 

43 

44 Returns 

45 ------- 

46 list of Orbits objs 

47 orbits associated with the list of input clusters 

48 

49 """ 

50 

51 logger.debug('Preparing input...') 

52 atoms = prepare_atoms(atom_list) 

53 clusters = prepare_clusters(cluster_list) 

54 rotations = prepare_rotations(rotation_matrices) 

55 translations = prepare_translations(translation_vectors) 

56 permutations = prepare_permutations(permutations) 

57 cell = prim.cell 

58 basis = prim.basis 

59 

60 logger.debug('Creating permutation map...') 

61 permutation_map, extended_atoms = \ 

62 get_permutation_map(atoms, rotations, translations, basis, symprec) 

63 

64 logger.debug('Creating orbits...') 

65 orbits = _get_orbits(permutation_map, extended_atoms, clusters, basis, 

66 cell, rotations, permutations) 

67 

68 return orbits 

69 

70 

71# All prepares are in case we changes some interface stuff 

72def prepare_permutations(permutations): 

73 perms = BiMap() 

74 for p in permutations: 

75 perms.append(tuple(p)) 

76 return perms 

77 

78 

79def prepare_atoms(atom_list): 

80 atoms = BiMap() 

81 for atom in atom_list: 

82 atoms.append(atom) 

83 return atoms 

84 

85 

86def prepare_clusters(cluster_list): 

87 clusters = BiMap() 

88 for cluster in cluster_list: 

89 clusters.append(tuple(cluster)) 

90 return clusters 

91 

92 

93def prepare_rotations(rotation_matrices): 

94 return rotation_matrices 

95 

96 

97def prepare_translations(translation_vectors): 

98 return translation_vectors 

99 

100 

101def get_permutation_map(atoms, rotations, translations, basis, symprec): 

102 

103 extended_atoms = atoms.copy() 

104 

105 permutation_map = np.zeros((len(atoms), len(rotations)), dtype=int) 

106 

107 scaled_positions = [atom.spos(basis) for atom in extended_atoms] 

108 

109 for sym_index, (R, T) in enumerate(zip(rotations, translations)): 

110 for atom_index, spos in enumerate(scaled_positions): 

111 

112 new_spos = np.dot(R, spos) + T 

113 new_atom = Atom.spos_to_atom(new_spos, basis, symprec) 

114 

115 if new_atom not in extended_atoms: 

116 extended_atoms.append(new_atom) 

117 

118 mapped_atom_index = extended_atoms.index(new_atom) 

119 permutation_map[atom_index, sym_index] = mapped_atom_index 

120 

121 return permutation_map, extended_atoms 

122 

123 

124# The internal function doing the outer loop over orbits 

125def _get_orbits(permutation_map, extended_atoms, clusters, 

126 basis, cell, 

127 rotations, permutations): 

128 cluster_is_found = [False] * len(clusters) 

129 orbits = [] 

130 for cluster_index, cluster in enumerate(clusters): 

131 if cluster_is_found[cluster_index]: 

132 continue 

133 

134 orbit = Orbit() 

135 

136 cluster_atoms = [extended_atoms[i] for i in cluster] 

137 

138 positions = [atom.pos(basis, cell) for atom in cluster_atoms] 

139 orbit.radius = get_geometrical_radius(positions) 

140 orbit.maximum_distance = get_maximum_distance(positions) 

141 orbit.order = len(cluster) 

142 

143 populate_orbit(orbit, permutations, clusters, 

144 cluster, 

145 permutation_map, extended_atoms, 

146 cluster_is_found) 

147 orbits.append(orbit) 

148 

149# bar.tick() 

150 return orbits 

151 

152 

153# Takes a cluster and generates all equivalent, translated clusters 

154def generate_translated_clusters(cluster, extended_atoms): 

155 transformed_cluster_atoms = [extended_atoms[i] for i in cluster] 

156 tested_offsets = set() 

157 for atom in transformed_cluster_atoms: 

158 offset = atom.offset 

159 if offset in tested_offsets: 

160 continue 

161 else: 

162 tested_offsets.add(offset) 

163 translated_cluster = [] 

164 for atom in transformed_cluster_atoms: 

165 new_offset = np.subtract(atom.offset, offset) 

166 new_atom = Atom(atom.site, new_offset) 

167 translated_cluster.append(extended_atoms.index(new_atom)) 

168 yield tuple(translated_cluster) 

169 

170 

171# Here is the actual categorization 

172def populate_orbit(orbit, permutations, clusters, 

173 cluster, 

174 permutation_map, extended_atoms, 

175 cluster_is_found): 

176 for sym_index in range(permutation_map.shape[1]): 

177 

178 of = OrientationFamily(sym_index) 

179 

180 transformed_cluster = permutation_map[cluster, sym_index] 

181 

182 for translated_cluster in generate_translated_clusters( 

183 transformed_cluster, extended_atoms): 

184 

185 argsort = tuple(np.argsort(translated_cluster)) 

186 perm_index = permutations.index(argsort) 

187 

188 translated_cluster = tuple(sorted(translated_cluster)) 

189 translated_cluster_index = clusters.index(translated_cluster) 

190 

191 if cluster == translated_cluster: 

192 if (sym_index, perm_index) not in orbit.eigensymmetries: 192 ↛ 195line 192 didn't jump to line 195, because the condition on line 192 was never false

193 orbit.eigensymmetries.append((sym_index, perm_index)) 

194 

195 if not cluster_is_found[translated_cluster_index]: 

196 cluster_is_found[translated_cluster_index] = True 

197 of.cluster_indices.append(translated_cluster_index) 

198 of.permutation_indices.append(perm_index) 

199 

200 if len(of.cluster_indices) > 0: 

201 orbit.orientation_families.append(of) 

202 

203 return orbit 

204 

205 

206class Orbit: 

207 """ 

208 This class serves as a container for storing data pertaining to an orbit. 

209 

210 Attributes 

211 ---------- 

212 orientation_families : list of OrientationFamily objs 

213 orientation families of the orbit 

214 eigensymmetries : list of tuples 

215 each eigensymmetry corresponds to a pair where the first index 

216 is the symmetry and the second is the permutation 

217 eigentensors : list(numpy.ndarray) 

218 decomposition of the force constant into symmetry elements 

219 """ 

220 

221 # TODO: Make properties of the parameters 

222 def __init__(self): 

223 self.orientation_families = [] 

224 self.eigensymmetries = [] 

225 self.eigentensors = [] 

226 self.radius = None 

227 self.order = None 

228 self.maximum_distance = None 

229 

230 @property 

231 def prototype_index(self): 

232 """int : index of cluster that serves as prototype for this orbit 

233 

234 In the code the first symmetry is always the identity so the first 

235 orientation family should always correspond to the prototype""" 

236 return self.orientation_families[0].cluster_indices[0] 

237 

238 def write(self, f): 

239 """Write a Orbit instance to a file. 

240 

241 Parameters 

242 ---------- 

243 f : str or file object 

244 name of input file (str) or stream to write to (file object) 

245 """ 

246 tar_file = tarfile.open(mode='w', fileobj=f) 

247 

248 # add all available attributes 

249 add_items_to_tarfile_pickle(tar_file, self.__dict__, 'attributes') 

250 

251 @staticmethod 

252 def read(f): 

253 """Load a ClusterSpace from file 

254 

255 Parameters 

256 ---------- 

257 f : string or file object 

258 name of input file (string) or stream to load from (file object) 

259 """ 

260 

261 orb = Orbit() 

262 tar_file = tarfile.open(mode='r', fileobj=f) 

263 

264 # read attributes pickle 

265 attributes = read_items_pickle(tar_file, 'attributes') 

266 for name, value in attributes.items(): 

267 orb.__setattr__(name, value) 

268 return orb 

269 

270 

271class OrientationFamily: 

272 """A container for storing information for a "family of orientations". 

273 

274 An orbit contains many clusters. Some of the clusters can be tranlsated 

275 onto each other and other must first be rotated. A set of clusters in the 

276 orbit which can all be translated onto each other are oriented in the same 

277 way and belongs to the same orientation family. The family is haracterized 

278 by the symmetry (rotation) which relates it to the prototype structure of 

279 the orbit. 

280 

281 Since the clusters are generally stored sorted the permutation information 

282 must also be stored. 

283 

284 Parameters 

285 ---------- 

286 symmetry_index : int 

287 The index of the symmetry corresponding to spglibs symmetry 

288 

289 Attributes 

290 ---------- 

291 symmetry_index : int 

292 The index of the symmetry corresponding to spglibs symmetry 

293 cluster_indices : list of ints 

294 The indices of the clusters belonging to this family 

295 permutation_indices : list of ints 

296 The indices of the permutation vector 

297 """ 

298 

299 def __init__(self, symmetry_index=None): 

300 self.symmetry_index = symmetry_index 

301 self.cluster_indices = [] 

302 self.permutation_indices = [] 

303 

304 def write(self, f): 

305 """ Write the object to file. 

306 

307 Parameters 

308 ---------- 

309 f : str or file object 

310 name of input file (str) or stream to write to (file object) 

311 """ 

312 pickle.dump(self, f) 

313 

314 @staticmethod 

315 def read(f): 

316 """Load a OrientationFamily object from a pickle file. 

317 

318 Parameters 

319 ---------- 

320 f : str or file object 

321 name of input file (str) or stream to load from (file object) 

322 

323 Returns 

324 ------- 

325 OrientationFamily object 

326 """ 

327 return pickle.load(f) 

328 

329 

330def get_geometrical_radius(positions): 

331 """Compute the geometrical size of a 3-dimensional point cloud. The 

332 geometrical size is defined as the average distance to the geometric 

333 center. 

334 

335 Parameters 

336 ---------- 

337 positions : list of 3-dimensional vectors 

338 positions of points in cloud 

339 

340 Returns 

341 ------- 

342 float 

343 geometric size of point cloud 

344 """ 

345 geometric_center = np.mean(positions, axis=0) 

346 return np.mean(np.sqrt(np.sum((positions - geometric_center)**2, axis=1))) 

347 

348 

349def get_maximum_distance(positions): 

350 """Compute the maximum distance between any two points in a 3-dimensional 

351 point cloud. This is equivalent to the "size" criterion used when imposing 

352 a certain (pair) cutoff criterion during construction of a set of clusters. 

353 

354 Parameters 

355 ---------- 

356 positions : list of 3-dimensional vectors 

357 positions of points in cloud 

358 

359 Returns 

360 ------- 

361 float 

362 maximum distance betwee any two points 

363 """ 

364 if len(positions) == 0: 364 ↛ 365line 364 didn't jump to line 365, because the condition on line 364 was never true

365 return 0.0 

366 max_distance = 0.0 

367 for pt1 in positions[:-1]: 

368 for pt2 in positions[1:]: 

369 max_distance = max(max_distance, np.linalg.norm(pt1 - pt2)) 

370 return max_distance