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

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

""" 

Contains the ClusterSpace object central to hiPhive 

""" 

 

import numpy as np 

import tarfile 

 

from ase.data import chemical_symbols 

from collections import OrderedDict 

from copy import deepcopy 

 

from .core.cluster_space_builder import build_cluster_space 

from .core.atoms import Atoms 

from .core.orbit import Orbit 

from .io.logging import logger 

from .io.read_write_files import (add_items_to_tarfile_pickle, 

add_items_to_tarfile_custom, 

add_list_to_tarfile_custom, 

read_items_pickle, 

read_list_custom) 

from .cluster_filter import Cutoffs, CutoffMaximumBody 

 

logger = logger.getChild('ClusterSpace') 

 

 

class ClusterSpace: 

# TODO: Inherit from base class to ease the use of FCP 

# TODO: Is this type of docstring correct? 

# TODO: Fix cutoffs parameter when new Cutoffs objects is implemented 

"""Primitive object handling cluster and force constants for a structure. 

 

Parameters 

---------- 

prototype_struture : ASE Atoms object 

prototype structure; spglib will be used to find a suitable cell based 

on this structure. 

cutoffs : list or Cutoffs obj 

cutoff radii for different orders starting with second order 

sum_rules : bool 

If True the aucostic sum rules will be enforced by constraining the 

parameters. 

symprec : float 

numerical precision that will be used for analyzing the symmetry (this 

parameter will be forwarded to 

`spglib <https://atztogo.github.io/spglib/>`_) 

length_scale : float 

This will be used as a normalization constant for the eigentensors 

 

Examples 

-------- 

To instantiate a ClusterSpace object one has to specify a prototype 

structure and cutoff radii for each cluster order that should be included. 

For example the following snippet will set up a ClusterSpace object for a 

BCC structure including second order terms up to a distance of 5 A and 

third order terms up to a distance of 4 A. 

 

>>> from ase.build import bulk 

>>> prim = bulk('W') 

>>> cs = ClusterSpace(prim, [5.0, 4.0]) 

""" 

# TODO: This class probably need some more documentation and a reference to 

# the thesis 

# TODO: Fix doc for n-body cutoff 

 

def __init__(self, prototype_structure, cutoffs, 

sum_rules=True, symprec=1e-5, length_scale=0.1): 

 

if not all(prototype_structure.pbc): 

raise ValueError('prototype_structure must have pbc.') 

 

if isinstance(cutoffs, Cutoffs): 

self._cutoffs = cutoffs 

if isinstance(cutoffs, list): 

self._cutoffs = CutoffMaximumBody(cutoffs, len(cutoffs) + 1) 

 

self._symprec = symprec 

self._length_scale = length_scale 

self._sum_rules = sum_rules 

 

self._atom_list = None 

self._cluster_list = None 

self._symmetry_dataset = None 

self._permutations = None 

self._prim = None 

self._orbits = None 

 

self._constraint_vectors = None 

# TODO: How to handle the constraint matrices? Should they even be 

# stored? 

self._constraint_matrices = None 

# Is this the best way or should the prim be instantiated separately? 

 

build_cluster_space(self, prototype_structure) 

 

# TODO: Should everything here be properties? deepcopy/ref etc.? 

# TODO: Docstrings for properties 

@property 

def number_of_dofs(self): 

"""int : number of free parameters in the model 

 

If the sum rules are not enforced the number of DOFs is the same as 

the total number of eigentensors in all orbits. 

""" 

return self._get_number_of_dofs() 

 

@property 

def cutoffs(self): 

""" Cutoffs object : cutoffs used for constructing the cluster space 

""" 

return deepcopy(self._cutoffs) 

 

@property 

def symprec(self): 

""" float : symprec value used when constructing the cluster space 

""" 

return self._symprec 

 

@property 

def sum_rules(self): 

""" bool : True if sum rules are enforced """ 

return self._sum_rules 

 

@property 

def length_scale(self): 

""" float : normalization constant of the force constants """ 

return self._length_scale 

 

@property 

def primitive_structure(self): 

""" ASE Atoms object : structure of the lattice """ 

return self._prim.copy() 

 

@property 

def spacegroup(self): 

""" str : space group of the lattice structure obtained from spglib """ 

return self._symmetry_dataset['international'] + ' ({})'.format( 

self._symmetry_dataset['number']) 

 

@property 

def wyckoff_sites(self): 

""" list : wyckoff sites in the primitive cell""" 

return self._symmetry_dataset['equivalent_atoms'] 

 

@property 

def rotation_matrices(self): 

""" list of 3x3 matrices : symmetry elements representing rotations """ 

return self._symmetry_dataset['rotations'].copy() 

 

@property 

def translation_vectors(self): 

""" list of 3-vectors : symmetry elements representin translations """ 

# TODO: bug incoming! 

return (self._symmetry_dataset['translations'] % 1).copy() 

 

@property 

def permutations(self): 

""" list of vectors : lookup for permutation references """ 

return deepcopy(self._permutations) 

 

@property 

def atom_list(self): 

""" BiMap object : atoms inside the cutoff relative to the of the 

center cell """ 

return self._atom_list 

 

@property 

def cluster_list(self): 

""" BiMap object : clusters possible within the cutoff """ 

return self._cluster_list 

 

@property 

def orbits(self): # TODO: add __getitem__ method 

""" list of Orbit objects : orbits associated with the lattice 

structure. """ 

return self._orbits 

 

def get_parameter_indices(self, order): 

""" 

Returns a list of the parameter indices associated with the requested 

order. 

 

Parameters 

---------- 

order : int 

order for which to return the parameter indices 

 

Returns 

------- 

list 

list of the parameter indices associated with the requested order 

 

Raises 

------ 

ValueError 

if the order is not included in the cluster space 

""" 

order = int(order) 

198 ↛ 199line 198 didn't jump to line 199, because the condition on line 198 was never true if order not in self.cutoffs.orders: 

raise ValueError('Order must be in {}'.format(self.cutoffs.orders)) 

min_param = 0 

max_param = 0 

for orbit in self.orbits: 

if orbit.order < order: 

min_param += len(orbit.eigentensors) 

max_param = min_param 

elif orbit.order == order: 

max_param += len(orbit.eigentensors) 

else: 

break 

rows, cols = self._cvs.nonzero() 

parameters = set() 

for r, c in zip(rows, cols): 

if min_param <= r < max_param: 

parameters.add(c) 

for r, c in zip(rows, cols): 

if c in parameters: 

assert min_param <= r < max_param, 'The internals are broken!' 

 

return sorted(parameters) 

 

@property 

def orbit_data(self): 

""" list : list of dictionaries containing detailed information for 

each orbit, e.g., cluster radius and atom types. 

""" 

data = [] 

p = 0 

for orbit_index, orbit in enumerate(self.orbits): 

d = {} 

d['index'] = orbit_index 

d['order'] = orbit.order 

d['radius'] = orbit.radius 

d['maximum_distance'] = orbit.maximum_distance 

d['number_of_clusters'] = len(orbit.orientation_families) 

d['eigentensors'] = orbit.eigentensors 

d['number_of_parameters'] = len(d['eigentensors']) 

 

types, wyckoff_sites = [], [] 

for atom_index in self.cluster_list[orbit.prototype_index]: 

atom = self.atom_list[atom_index] 

types.append(self.primitive_structure.numbers[atom.site]) 

wyckoff_sites.append(self.wyckoff_sites[atom.site]) 

d['prototype_cluster'] = self.cluster_list[orbit.prototype_index] 

d['prototype_atom_types'] = tuple(types) 

d['prototype_wyckoff_sites'] = tuple(wyckoff_sites) 

 

d['geometrical_order'] = len(set(d['prototype_cluster'])) 

d['parameter_indices'] = np.arange(p, p + len(orbit.eigentensors)) 

 

p += len(orbit.eigentensors) 

data.append(d) 

 

return data 

 

def print_orbits(self): 

""" Prints a list of all orbits. """ 

orbits = self.orbit_data 

 

def str_orbit(index, orbit): 

elements = ' '.join(chemical_symbols[n] for n in 

orbit['prototype_atom_types']) 

fields = OrderedDict([ 

('index', '{:^3}'.format(index)), 

('order', '{:^3}'.format(orbit['order'])), 

('elements', '{:^18}'.format(elements)), 

('radius', '{:^8.4f}'.format(orbit['radius'])), 

('prototype', '{:^18}' 

.format(str(orbit['prototype_cluster']))), 

('clusters', '{:^4}'.format(orbit['number_of_clusters'])), 

('parameters', '{:^3}'.format(len(orbit['eigentensors']))), 

]) 

 

s = [] 

for name, value in fields.items(): 

n = max(len(name), len(value)) 

if index < 0: 

s += ['{s:^{n}}'.format(s=name, n=n)] 

else: 

s += ['{s:^{n}}'.format(s=value, n=n)] 

return ' | '.join(s) 

 

# table header 

width = max(len(str_orbit(-1, orbits[-1])), 

len(str_orbit(0, orbits[-1]))) 

print(' List of Orbits '.center(width, '=')) 

print(str_orbit(-1, orbits[0])) 

print(''.center(width, '-')) 

 

# table body 

for i, orbit in enumerate(orbits): 

print(str_orbit(i, orbit)) 

print(''.center(width, '=')) 

 

def _get_number_of_dofs(self): 

""" Returns the number of degrees of freedom. """ 

return self._cvs.shape[1] 

 

def __str__(self): 

 

def str_order(order, header=False): 

formats = {'order': '{:2}', 

'num_orbits': '{:5}', 

'num_params': '{:5}', 

'num_clusters': '{:5}'} 

s = [] 

for name, value in order.items(): 

str_repr = formats[name].format(value) 

n = max(len(name), len(str_repr)) 

if header: 

s += ['{s:^{n}}'.format(s=name, n=n)] 

else: 

s += ['{s:^{n}}'.format(s=str_repr, n=n)] 

return ' | '.join(s) 

 

# collect data 

orbits = self.orbit_data 

orders = OrderedDict( 

(order, OrderedDict(order=order)) for order in self.cutoffs.orders) 

for orbit in orbits: 

k = orbit['order'] 

orders[k]['num_orbits'] = orders[k].get('num_orbits', 0) + 1 

orders[k]['num_clusters'] = orders[k].get('num_clusters', 0) + \ 

orbit['number_of_clusters'] 

 

# basic information 

prototype = orders[list(orders.keys())[-1]] 

n = max(len(str_order(prototype)), 54) 

s = [] 

s.append(' Cluster Space '.center(n, '=')) 

data = [('Spacegroup', self.spacegroup), 

('symprec', self.symprec), 

('Sum rules', self.sum_rules), 

('Length scale', self.length_scale), 

('Cutoffs', self.cutoffs), 

('Cell', self.primitive_structure.cell), 

('Basis', self.primitive_structure.basis), 

('Numbers', self.primitive_structure.numbers), 

('Total number of orbits', len(orbits)), 

('Total number of clusters', 

sum([order['num_clusters'] for order in orders.values()])), 

('Total number of parameters', 

self._get_number_of_dofs() 

)] 

for field, value in data: 

if str(value).count('\n') > 1: 

s.append('{:26} :\n{}'.format(field, value)) 

else: 

s.append('{:26} : {}'.format(field, value)) 

 

# table header 

s.append(''.center(n, '-')) 

s.append(str_order(prototype, header=True)) 

s.append(''.center(n, '-')) 

for order in orders.values(): 

s.append(str_order(order).rstrip()) 

s.append(''.center(n, '=')) 

return '\n'.join(s) 

 

def __repr__(self): 

s = 'ClusterSpace({!r}, {!r}, {!r}, {!r}, {!r})' 

return s.format(self.primitive_structure, self.cutoffs, self.sum_rules, 

self.symprec, self.length_scale) 

 

def _map_parameters(self, parameters): 

""" Maps irreducible parameters to the real parameters associated with 

the eigentensors. 

""" 

368 ↛ 369line 368 didn't jump to line 369, because the condition on line 368 was never true if len(parameters) != self.number_of_dofs: 

raise ValueError('Invalid number of parameters, please provide {} ' 

'parameters'.format(self.number_of_dofs)) 

return self._cvs.dot(parameters) 

 

def copy(self): 

return deepcopy(self) 

 

def write(self, fileobj): 

""" Saves to file or to a file-like object. 

 

The instance is saved into a custom format based on tar-files. The 

resulting file will be a valid tar file and can be browsed by by a tar 

reader. The included objects are themself either pickles, npz or other 

tars. 

 

Parameters 

---------- 

fileobj : str or file-like object 

If the input is a string a tar archive will be created in the 

current directory. Otherwise the input must be a valid file 

like object. 

""" 

# Create a tar archive 

if isinstance(fileobj, str): 

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

else: 

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

 

# Attributes in pickle format 

pickle_attributes = ['_symprec', '_length_scale', '_sum_rules', 

'_symmetry_dataset', '_permutations', 

'_atom_list', '_cluster_list'] 

items_pickle = dict() 

for attribute in pickle_attributes: 

items_pickle[attribute] = self.__getattribute__(attribute) 

add_items_to_tarfile_pickle(tar_file, items_pickle, 'attributes') 

 

# Constraint matrices and vectors in hdf5 format 

items = dict(cvs=self._cvs) 

add_items_to_tarfile_pickle(tar_file, items, 'constraint_vectors') 

 

# Cutoffs and prim with their builtin write/read functions 

items_custom = {'_cutoffs': self._cutoffs, '_prim': self._prim} 

add_items_to_tarfile_custom(tar_file, items_custom) 

 

# Orbits 

add_list_to_tarfile_custom(tar_file, self._orbits, 'orbits') 

add_list_to_tarfile_custom(tar_file, self._dropped_orbits, 

'dropped_orbits') 

 

# Done! 

tar_file.close() 

 

def read(f): 

""" Reads a ClusterSpace instance from file. 

 

Parameters 

---------- 

f : string or file object 

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

""" 

 

# Instantiate empty cs obj. 

cs = ClusterSpace.__new__(ClusterSpace) 

 

# Load from file on disk or file-like 

if type(f) is str: 

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

else: 

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

 

# Attributes 

attributes = read_items_pickle(tar_file, 'attributes') 

for name, value in attributes.items(): 

cs.__setattr__(name, value) 

 

# Load the constraint matrices into their dict 

items = read_items_pickle(tar_file, 'constraint_vectors') 

cs._cvs = items['cvs'] 

 

# Cutoffs and prim via custom save funcs 

fileobj = tar_file.extractfile('_cutoffs') 

cs._cutoffs = Cutoffs.read(fileobj) 

 

fileobj = tar_file.extractfile('_prim') 

cs._prim = Atoms.read(fileobj) 

 

# Orbits are stored in a separate archive 

cs._orbits = read_list_custom(tar_file, 'orbits', Orbit.read) 

cs._dropped_orbits = read_list_custom( 

tar_file, 'dropped_orbits', Orbit.read) 

 

return cs