Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

""" 

Contains the Orbit class which hold onformation about equivalent clusters. 

""" 

import pickle 

import numpy as np 

import tarfile 

 

from .utilities import BiMap 

from .atoms import Atom 

from ..input_output.logging_tools import logger 

from ..input_output.read_write_files import (add_items_to_tarfile_hdf5, 

add_items_to_tarfile_pickle, 

add_list_to_tarfile_custom, 

read_items_hdf5, 

read_items_pickle, 

read_list_custom) 

logger = logger.getChild('orbits') 

 

 

# This is the interface accessible for cluster_space 

def get_orbits(cluster_list, atom_list, rotation_matrices, translation_vectors, 

permutations, prim, symprec): 

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

configuration. 

 

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

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

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

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

respective set of equivalent clusters. 

 

Parameters 

---------- 

cluster_list : BiMap object 

a list of clusters 

atom_list : BiMap object 

a list of atoms in a supercell 

rotation_matrices : list of NumPy (3,3) arrays 

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

translation_vectors : list of NumPy (3) arrays 

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

permutations : list of permutations 

lookup table for permutations 

prim : hiPhive Atoms object 

primitive structure 

 

Returns 

------- 

list of Orbits objs 

orbits associated with the list of input clusters 

 

""" 

 

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

atoms = prepare_atoms(atom_list) 

clusters = prepare_clusters(cluster_list) 

rotations = prepare_rotations(rotation_matrices) 

translations = prepare_translations(translation_vectors) 

permutations = prepare_permutations(permutations) 

cell = prim.cell 

basis = prim.basis 

 

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

permutation_map, extended_atoms = \ 

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

 

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

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

cell, rotations, permutations) 

 

return orbits 

 

 

# All prepares are in case we changes some interface stuff 

def prepare_permutations(permutations): 

perms = BiMap() 

for p in permutations: 

perms.append(tuple(p)) 

return perms 

 

 

def prepare_atoms(atom_list): 

atoms = BiMap() 

for atom in atom_list: 

atoms.append(atom) 

return atoms 

 

 

def prepare_clusters(cluster_list): 

clusters = BiMap() 

for cluster in cluster_list: 

clusters.append(tuple(cluster)) 

return clusters 

 

 

def prepare_rotations(rotation_matrices): 

return rotation_matrices 

 

 

def prepare_translations(translation_vectors): 

return translation_vectors 

 

 

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

 

extended_atoms = atoms.copy() 

 

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

 

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

 

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

for atom_index, spos in enumerate(scaled_positions): 

 

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

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

 

if new_atom not in extended_atoms: 

extended_atoms.append(new_atom) 

 

mapped_atom_index = extended_atoms.index(new_atom) 

permutation_map[atom_index, sym_index] = mapped_atom_index 

 

return permutation_map, extended_atoms 

 

 

# The internal function doing the outer loop over orbits 

def _get_orbits(permutation_map, extended_atoms, clusters, 

basis, cell, 

rotations, permutations): 

cluster_is_found = [False] * len(clusters) 

orbits = [] 

for cluster_index, cluster in enumerate(clusters): 

if cluster_is_found[cluster_index]: 

continue 

 

orbit = Orbit() 

 

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

 

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

orbit.radius = get_geometrical_radius(positions) 

orbit.maximum_distance = get_maximum_distance(positions) 

orbit.order = len(cluster) 

 

populate_orbit(orbit, permutations, clusters, 

cluster, 

permutation_map, extended_atoms, 

cluster_is_found) 

orbits.append(orbit) 

 

# bar.tick() 

return orbits 

 

 

# Takes a cluster and generates all equivalent, translated clusters 

def generate_translated_clusters(cluster, extended_atoms): 

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

tested_offsets = set() 

for atom in transformed_cluster_atoms: 

offset = atom.offset 

if offset in tested_offsets: 

continue 

else: 

tested_offsets.add(offset) 

translated_cluster = [] 

for atom in transformed_cluster_atoms: 

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

new_atom = Atom(atom.site, new_offset) 

translated_cluster.append(extended_atoms.index(new_atom)) 

yield tuple(translated_cluster) 

 

 

# Here is the actual categorization 

def populate_orbit(orbit, permutations, clusters, 

cluster, 

permutation_map, extended_atoms, 

cluster_is_found): 

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

 

of = OrientationFamily(sym_index) 

 

transformed_cluster = permutation_map[cluster, sym_index] 

 

for translated_cluster in generate_translated_clusters( 

transformed_cluster, extended_atoms): 

 

argsort = tuple(np.argsort(translated_cluster)) 

perm_index = permutations.index(argsort) 

 

translated_cluster = tuple(sorted(translated_cluster)) 

translated_cluster_index = clusters.index(translated_cluster) 

 

if cluster == translated_cluster: 

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

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

 

if not cluster_is_found[translated_cluster_index]: 

cluster_is_found[translated_cluster_index] = True 

of.cluster_indices.append(translated_cluster_index) 

of.permutation_indices.append(perm_index) 

 

if len(of.cluster_indices) > 0: 

orbit.orientation_families.append(of) 

 

return orbit 

 

 

class Orbit: 

""" 

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

 

Attributes 

---------- 

orientation_families : list of OrientationFamily objs 

orientation families of the orbit 

eigensymmetries : list of tuples 

each eigensymmetry corresponds to a pair where the first index 

is the symmetry and the second is the permutation 

eigentensors : list(numpy.ndarray) 

decomposition of the force constant into symmetry elements 

""" 

 

# TODO: Make properties of the parameters 

def __init__(self): 

self.orientation_families = [] 

self.eigensymmetries = [] 

self.eigentensors = [] 

 

@property 

def prototype_index(self): 

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

 

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

orientation family should always correspond to the prototype""" 

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

 

def write(self, f): 

"""Write a Orbit instance to a file. 

 

Parameters 

---------- 

f : str or file object 

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

""" 

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

 

# add eigentensors as NumPy array 

items_hdf5 = dict(eigentensors=np.array(self.eigentensors), 

) 

add_items_to_tarfile_hdf5(tar_file, items_hdf5, 'eigentensors') 

 

# add eigensymmetries as list 

items_pickle = dict(eigensymmetries=self.eigensymmetries, 

order=self.order, radius=self.radius, 

maximum_distance=self.maximum_distance) 

add_items_to_tarfile_pickle(tar_file, items_pickle, 'attributes') 

 

# add orientation families 

add_list_to_tarfile_custom(tar_file, self.orientation_families, 

'orientation_families') 

tar_file.close() 

 

@staticmethod 

def read(f): 

"""Load a ClusterSpace from file 

 

Parameters 

---------- 

f : string or file object 

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

""" 

 

orb = Orbit() 

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

 

# read eigentensors hdf5 

items_hdf5 = read_items_hdf5(tar_file, 'eigentensors') 

orb.eigentensors = [et for et in items_hdf5['eigentensors']] 

 

# read attributes pickle 

attributes = read_items_pickle(tar_file, 'attributes') 

for name, value in attributes.items(): 

orb.__setattr__(name, value) 

 

# read orientation families 

ofs = read_list_custom(tar_file, 'orientation_families', 

OrientationFamily.read) 

orb.orientation_families = ofs 

return orb 

 

 

class OrientationFamily: 

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

 

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

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

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

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

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

the orbit. 

 

Since the clusters are generally stored sorted the permutation information 

must also be stored. 

 

Parameters 

---------- 

symmetry_index : int 

The index of the symmetry corresponding to spglibs symmetry 

 

Attributes 

---------- 

symmetry_index : int 

The index of the symmetry corresponding to spglibs symmetry 

cluster_indices : list of ints 

The indices of the clusters belonging to this family 

permutation_indices : list of ints 

The indices of the permutation vector 

""" 

 

def __init__(self, symmetry_index=None): 

self.symmetry_index = symmetry_index 

self.cluster_indices = [] 

self.permutation_indices = [] 

 

def write(self, f): 

""" Write the object to file. 

 

Parameters 

---------- 

f : str or file object 

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

""" 

pickle.dump(self, f) 

 

@staticmethod 

def read(f): 

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

 

Parameters 

---------- 

f : str or file object 

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

 

Returns 

------- 

OrientationFamily object 

""" 

return pickle.load(f) 

 

 

def get_geometrical_radius(positions): 

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

geometrical size is defined as the average distance to the geometric 

center. 

 

Parameters 

---------- 

positions : list of 3-dimensional vectors 

positions of points in cloud 

 

Returns 

------- 

float 

geometric size of point cloud 

""" 

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

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

 

 

def get_maximum_distance(positions): 

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

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

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

 

Parameters 

---------- 

positions : list of 3-dimensional vectors 

positions of points in cloud 

 

Returns 

------- 

float 

maximum distance betwee any two points 

""" 

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

return 0.0 

max_distance = 0.0 

for pt1 in positions[:-1]: 

for pt2 in positions[1:]: 

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

return max_distance