Coverage for hiphive/cluster

1"""

2Contains the ClusterSpace object central to hiPhive

3"""

5import numpy as np

6import tarfile

8from ase.data import chemical_symbols

9from collections import OrderedDict

10from copy import deepcopy

12from .core.cluster_space_builder import build_cluster_space

13from .core.atoms import Atoms

14from .core.orbits import Orbit

15from .cluster_space_data import ClusterSpaceData

16from .input_output.logging_tools import logger

17from .input_output.read_write_files import (add_items_to_tarfile_pickle,

18 add_items_to_tarfile_custom,

19 add_list_to_tarfile_custom,

20 read_items_pickle,

21 read_list_custom)

22from .cutoffs import Cutoffs, CutoffMaximumBody, BaseClusterFilter

24from .config import Config

25logger = logger.getChild('ClusterSpace')

28class ClusterSpace:

29 """Primitive object for handling clusters and force constants of a structure.

31 Parameters

32 ----------

33 prototype_structure : ase.Atoms

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

35 on this structure unless the cell is already a primitive cell.

36 cutoffs : list or Cutoffs

37 cutoff radii for different orders starting with second order

38 cluster_filter : ClusterFilter

39 accepts a subclass of hiphive.filters.BaseClusterFilter to further

40 control which orbits to include.

41 config : Config object

42 a configuration object that holds information on how the cluster space

43 should be built, e.g., values for tolerances and specifications

44 regarding the handling of acoustic sum rules; if ``config`` is

45 not given then the keyword arguments that follow below can be

46 used for configuration.

47 acoustic_sum_rules : bool

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

49 parameters.

50 symprec : float

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

52 parameter will be forwarded to

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

54 length_scale : float

55 This will be used as a normalization constant for the eigentensors

57 Examples

58 --------

60 To instantiate a :class:`ClusterSpace` object one has to specify a

61 prototype structure and cutoff radii for each cluster order that

62 should be included. For example the following snippet will set up

63 a :class:`ClusterSpace` object for a body-centered-cubic (BCC)

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

65 third order terms up to a distance of 4 A.

67 >>> from ase.build import bulk

68 >>> from hiphive import ClusterSpace

69 >>> prim = bulk('W')

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

72 """

73 # TODO: This class probably needs some more documentation

75 def __init__(self, prototype_structure, cutoffs, config=None,

76 cluster_filter=None, **kwargs):

78 if not all(prototype_structure.pbc):

79 raise ValueError('prototype_structure must have pbc.')

81 if isinstance(cutoffs, Cutoffs):

82 self._cutoffs = cutoffs

83 elif isinstance(cutoffs, list):

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

85 else:

86 raise TypeError('cutoffs must be a list or a Cutoffs object')

88 if config is None:

89 config = Config(**kwargs)

90 else:

91 if not isinstance(config, Config): 91 ↛ 92line 91 didn't jump to line 92, because the condition on line 91 was never true

92 raise TypeError('config kw must be of type {}'.format(Config))

93 if len(kwargs) > 0:

94 raise ValueError('use either Config or kwargs, not both')

95 self._config = config

97 if cluster_filter is None:

98 self._cluster_filter = BaseClusterFilter()

99 else:

100 self._cluster_filter = cluster_filter

101

102 self._atom_list = None

103 self._cluster_list = None

104 self._symmetry_dataset = None

105 self._permutations = None

106 self._prim = None

107 self._orbits = None

108

109 self._constraint_vectors = None

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

111 # stored?

112 self._constraint_matrices = None

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

114

115 build_cluster_space(self, prototype_structure)

116 self.summary = ClusterSpaceData(self)

117

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

119 @property

120 def n_dofs(self):

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

122

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

124 the total number of eigentensors in all orbits.

125 """

126 return self._get_n_dofs()

127

128 @property

129 def cutoffs(self):

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

131 return deepcopy(self._cutoffs)

132

133 @property

134 def symprec(self):

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

136 return self._config['symprec']

137

138 @property

139 def acoustic_sum_rules(self):

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

141 return self._config['acoustic_sum_rules']

142

143 @property

144 def length_scale(self):

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

146 return self._config['length_scale']

147

148 @property

149 def primitive_structure(self):

150 """ ase.Atoms : structure of the lattice """

151 return self._prim.copy()

152

153 @property

154 def spacegroup(self):

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

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

157 self._symmetry_dataset['number'])

158

159 @property

160 def wyckoff_sites(self):

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

162 return self._symmetry_dataset['equivalent_atoms']

163

164 @property

165 def rotation_matrices(self):

166 """ list(numpy.ndarray) : symmetry elements (`3x3` matrices)

167 representing rotations """

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

169

170 @property

171 def translation_vectors(self):

172 """ list(numpy.ndarray) : symmetry elements representing

173 translations """

174 # TODO: bug incoming!

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

176

177 @property

178 def permutations(self):

179 """ list(numpy.ndarray) : lookup for permutation references """

180 return deepcopy(self._permutations)

181

182 @property

183 def atom_list(self):

184 """ BiMap : atoms inside the cutoff relative to the of the center cell

185 """

186 return self._atom_list

187

188 @property

189 def cluster_list(self):

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

191 return self._cluster_list

192

193 @property

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

195 """ list(Orbit) : orbits associated with the lattice structure. """

196 return self._orbits

197

198 @property

199 def orbit_data(self):

200 """ list(dict) : detailed information for each orbit, e.g., cluster

201 radius and atom types.

202 """

203 data = []

204 p = 0

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

206 d = {}

207 d['index'] = orbit_index

208 d['order'] = orbit.order

209 d['radius'] = orbit.radius

210 d['maximum_distance'] = orbit.maximum_distance

211 d['n_clusters'] = len(orbit.orientation_families)

212 d['eigentensors'] = orbit.eigentensors

213 d['n_parameters'] = len(d['eigentensors'])

214

215 types, wyckoff_sites = [], []

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

217 atom = self.atom_list[atom_index]

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

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

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

221 d['prototype_atom_types'] = tuple(types)

222 d['prototype_wyckoff_sites'] = tuple(wyckoff_sites)

223

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

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

226

227 p += len(orbit.eigentensors)

228 data.append(d)

229

230 return data

231

232 def get_parameter_indices(self, order):

233 """

234 Returns a list of the parameter indices associated with the requested

235 order.

236

237 Parameters

238 ----------

239 order : int

240 order for which to return the parameter indices

241

242 Returns

243 -------

244 list(int)

245 list of parameter indices associated with the requested order

246

247 Raises

248 ------

249 ValueError

250 if the order is not included in the cluster space

251 """

252 order = int(order)

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

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

255 min_param = 0

256 max_param = 0

257 for orbit in self.orbits:

258 if orbit.order < order:

259 min_param += len(orbit.eigentensors)

260 max_param = min_param

261 elif orbit.order == order:

262 max_param += len(orbit.eigentensors)

263 else:

264 break

265 rows, cols = self._cvs.nonzero()

266 parameters = set()

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

268 if min_param <= r < max_param:

269 parameters.add(c)

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

271 if c in parameters:

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

273

274 return sorted(parameters)

275

276 def get_n_dofs_by_order(self, order):

277 """ Returns number of degrees of freedom for the given order.

278

279 Parameters

280 ----------

281 order : int

282 order for which to return the number of dofs

283

284 Returns

285 -------

286 int

287 number of degrees of freedom

288 """

289 return len(self.get_parameter_indices(order=order))

290

291 def _get_n_dofs(self):

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

293 return self._cvs.shape[1]

294

295 def _map_parameters(self, parameters):

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

297 the eigentensors.

298 """

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

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

301 'parameters'.format(self.n_dofs))

302 return self._cvs.dot(parameters)

303

304 def print_tables(self):

305 """ Prints information concerning the underlying cluster space to stdout,including,

306 e.g., the number of cluster, orbits, and parameters by order and number of bodies. """

307 self.summary.print_tables()

308

309 def print_orbits(self):

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

311 orbits = self.orbit_data

312

313 def str_orbit(index, orbit):

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

315 orbit['prototype_atom_types'])

316 fields = OrderedDict([

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

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

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

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

321 ('prototype', '{:^18}'.format(str(orbit['prototype_cluster']))),

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

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

324 ])

325

326 s = []

327 for name, value in fields.items():

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

329 if index < 0:

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

331 else:

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

333 return ' | '.join(s)

334

335 # table header

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

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

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

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

340

341 # table body

342 for i, orbit in enumerate(orbits):

343 print(str_orbit(i, orbit))

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

345

346 def __str__(self):

347

348 def str_order(order, header=False):

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

350 'n_orbits': '{:5}',

351 'n_clusters': '{:5}'}

352 s = []

353 for name, value in order.items():

354 str_repr = formats[name].format(value)

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

356 if header:

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

358 else:

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

360 return ' | '.join(s)

361

362 # collect data

363 orbits = self.orbit_data

364 orders = self.cutoffs.orders

365

366 order_data = {o: dict(order=o, n_orbits=0, n_clusters=0) for o in orders}

367 for orbit in orbits:

368 o = orbit['order']

369 order_data[o]['n_orbits'] += 1

370 order_data[o]['n_clusters'] += orbit['n_clusters']

371

372 # prototype with max order to find column width

373 max_order = max(orders)

374 prototype = order_data[max_order]

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

376

377 # basic information

378 s = []

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

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

381 ('symprec', self.symprec),

382 ('Sum rules', self.acoustic_sum_rules),

383 ('Length scale', self.length_scale),

384 ('Cutoffs', self.cutoffs),

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

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

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

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

389 ('Total number of clusters',

390 sum([order_data[order]['n_clusters'] for order in orders])),

391 ('Total number of parameters', self._get_n_dofs()

392 )]

393 for field, value in data:

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

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

396 else:

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

398

399 # table header

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

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

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

403 for order in orders:

404 s.append(str_order(order_data[order]).rstrip())

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

406 return '\n'.join(s)

407

408 def __repr__(self):

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

410 return s.format(self.primitive_structure, self.cutoffs, self._config)

411

412 def copy(self):

413 return deepcopy(self)

414

415 def write(self, fileobj):

416 """ Writes cluster space to file.

417

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

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

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

421 tars.

422

423 Parameters

424 ----------

425 fileobj : str or file-like object

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

427 current directory. Otherwise the input must be a valid file

428 like object.

429 """

430 # Create a tar archive

431 if isinstance(fileobj, str):

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

433 else:

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

435

436 # Attributes in pickle format

437 pickle_attributes = ['_config',

438 '_symmetry_dataset', '_permutations',

439 '_atom_list', '_cluster_list']

440 items_pickle = dict()

441 for attribute in pickle_attributes:

442 items_pickle[attribute] = self.__getattribute__(attribute)

443 add_items_to_tarfile_pickle(tar_file, items_pickle, 'attributes')

444

445 # Constraint matrices and vectors in hdf5 format

446 items = dict(cvs=self._cvs)

447 add_items_to_tarfile_pickle(tar_file, items, 'constraint_vectors')

448

449 # Cutoffs and prim with their builtin write/read functions

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

451 add_items_to_tarfile_custom(tar_file, items_custom)

452

453 # Orbits

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

455 add_list_to_tarfile_custom(tar_file, self._dropped_orbits,

456 'dropped_orbits')

457

458 # Done!

459 tar_file.close()

460

461 def read(f):

462 """ Reads a cluster space from file.

463

464 Parameters

465 ----------

466 f : str or file object

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

468 """

469

470 # Instantiate empty cs obj.

471 cs = ClusterSpace.__new__(ClusterSpace)

472

473 # Load from file on disk or file-like

474 if type(f) is str:

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

476 else:

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

478

479 # Attributes

480 attributes = read_items_pickle(tar_file, 'attributes')

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

482 cs.__setattr__(name, value)

483

484 # Load the constraint matrices into their dict

485 items = read_items_pickle(tar_file, 'constraint_vectors')

486 cs._cvs = items['cvs']

487

488 # Cutoffs and prim via custom save funcs

489 fileobj = tar_file.extractfile('_cutoffs')

490 cs._cutoffs = Cutoffs.read(fileobj)

491

492 fileobj = tar_file.extractfile('_prim')

493 cs._prim = Atoms.read(fileobj)

494

495 # Orbits are stored in a separate archive

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

497 cs._dropped_orbits = read_list_custom(

498 tar_file, 'dropped_orbits', Orbit.read)

499

500 tar_file.close()

501

502 # create summary object based on CS

503 cs.summary = ClusterSpaceData(cs)

504

505 # Done!

506 return cs

Coverage for hiphive/cluster_space.py: 98%

234 statements