Coverage for /builds/kinetik161/ase/ase/spacegroup/spacegroup.py: 92.17%

1# Copyright (C) 2010, Jesper Friis 

2# (see accompanying license files for details). 

3"""Definition of the Spacegroup class. 

4 

5This module only depends on NumPy and the space group database. 

6""" 

7 

8import os 

9import warnings 

10from functools import total_ordering 

11from typing import Union 

12 

13import numpy as np 

14 

15__all__ = ['Spacegroup'] 

16 

17 

18class SpacegroupError(Exception): 

19 """Base exception for the spacegroup module.""" 

20 

21 

22class SpacegroupNotFoundError(SpacegroupError): 

23 """Raised when given space group cannot be found in data base.""" 

24 

25 

26class SpacegroupValueError(SpacegroupError): 

27 """Raised when arguments have invalid value.""" 

28 

29 

30# Type alias 

31_SPACEGROUP = Union[int, str, 'Spacegroup'] 

32 

33 

34@total_ordering 

35class Spacegroup: 

36 """A space group class. 

37 

38 The instances of Spacegroup describes the symmetry operations for 

39 the given space group. 

40 

41 Example: 

42 

43 >>> from ase.spacegroup import Spacegroup 

44 >>> 

45 >>> sg = Spacegroup(225) 

46 >>> print('Space group', sg.no, sg.symbol) 

47 Space group 225 F m -3 m 

48 >>> sg.scaled_primitive_cell 

49 array([[ 0. , 0.5, 0.5], 

50 [ 0.5, 0. , 0.5], 

51 [ 0.5, 0.5, 0. ]]) 

52 >>> sites, kinds = sg.equivalent_sites([[0,0,0]]) 

53 >>> sites 

54 array([[ 0. , 0. , 0. ], 

55 [ 0. , 0.5, 0.5], 

56 [ 0.5, 0. , 0.5], 

57 [ 0.5, 0.5, 0. ]]) 

58 """ 

59 no = property( 

60 lambda self: self._no, 

61 doc='Space group number in International Tables of Crystallography.') 

62 symbol = property( 

63 lambda self: self._symbol, 

64 doc='Hermann-Mauguin (or international) symbol for the space group.') 

65 setting = property(lambda self: self._setting, 

66 doc='Space group setting. Either one or two.') 

67 lattice = property(lambda self: self._symbol[0], 

68 doc="""Lattice type: 

69 

70 P primitive 

71 I body centering, h+k+l=2n 

72 F face centering, h,k,l all odd or even 

73 A,B,C single face centering, k+l=2n, h+l=2n, h+k=2n 

74 R rhombohedral centering, -h+k+l=3n (obverse); h-k+l=3n (reverse) 

75 """) 

76 centrosymmetric = property(lambda self: self._centrosymmetric, 

77 doc='Whether a center of symmetry exists.') 

78 scaled_primitive_cell = property( 

79 lambda self: self._scaled_primitive_cell, 

80 doc='Primitive cell in scaled coordinates as a matrix with the ' 

81 'primitive vectors along the rows.') 

82 reciprocal_cell = property( 

83 lambda self: self._reciprocal_cell, 

84 doc='Tree Miller indices that span all kinematically non-forbidden ' 

85 'reflections as a matrix with the Miller indices along the rows.') 

86 nsubtrans = property(lambda self: len(self._subtrans), 

87 doc='Number of cell-subtranslation vectors.') 

88 

89 def _get_nsymop(self): 

90 """Returns total number of symmetry operations.""" 

91 if self.centrosymmetric: 

92 return 2 * len(self._rotations) * len(self._subtrans) 

93 else: 

94 return len(self._rotations) * len(self._subtrans) 

95 

96 nsymop = property(_get_nsymop, doc='Total number of symmetry operations.') 

97 subtrans = property( 

98 lambda self: self._subtrans, 

99 doc='Translations vectors belonging to cell-sub-translations.') 

100 rotations = property( 

101 lambda self: self._rotations, 

102 doc='Symmetry rotation matrices. The invertions are not included ' 

103 'for centrosymmetrical crystals.') 

104 translations = property( 

105 lambda self: self._translations, 

106 doc='Symmetry translations. The invertions are not included ' 

107 'for centrosymmetrical crystals.') 

108 

109 def __init__(self, spacegroup: _SPACEGROUP, setting=1, datafile=None): 

110 """Returns a new Spacegroup instance. 

111 

112 Parameters: 

113 

114 spacegroup : int | string | Spacegroup instance 

115 The space group number in International Tables of 

116 Crystallography or its Hermann-Mauguin symbol. E.g. 

117 spacegroup=225 and spacegroup='F m -3 m' are equivalent. 

118 setting : 1 | 2 

119 Some space groups have more than one setting. `setting` 

120 determines Which of these should be used. 

121 datafile : None | string 

122 Path to database file. If `None`, the the default database 

123 will be used. 

124 """ 

125 if isinstance(spacegroup, Spacegroup): 

126 for k, v in spacegroup.__dict__.items(): 

127 setattr(self, k, v) 

128 return 

129 if not datafile: 

130 datafile = get_datafile() 

131 with open(datafile) as fd: 

132 _read_datafile(self, spacegroup, setting, fd) 

133 

134 def __repr__(self): 

135 return 'Spacegroup(%d, setting=%d)' % (self.no, self.setting) 

136 

137 def todict(self): 

138 return {'number': self.no, 'setting': self.setting} 

139 

140 def __str__(self): 

141 """Return a string representation of the space group data in 

142 the same format as found the database.""" 

143 retval = [] 

144 # no, symbol 

145 retval.append('%-3d %s\n' % (self.no, self.symbol)) 

146 # setting 

147 retval.append(' setting %d\n' % (self.setting)) 

148 # centrosymmetric 

149 retval.append(' centrosymmetric %d\n' % (self.centrosymmetric)) 

150 # primitive vectors 

151 retval.append(' primitive vectors\n') 

152 for i in range(3): 

153 retval.append(' ') 

154 for j in range(3): 

155 retval.append(' %13.10f' % (self.scaled_primitive_cell[i, j])) 

156 retval.append('\n') 

157 # primitive reciprocal vectors 

158 retval.append(' reciprocal vectors\n') 

159 for i in range(3): 

160 retval.append(' ') 

161 for j in range(3): 

162 retval.append(' %3d' % (self.reciprocal_cell[i, j])) 

163 retval.append('\n') 

164 # sublattice 

165 retval.append(' %d subtranslations\n' % self.nsubtrans) 

166 for i in range(self.nsubtrans): 

167 retval.append(' ') 

168 for j in range(3): 

169 retval.append(' %13.10f' % (self.subtrans[i, j])) 

170 retval.append('\n') 

171 # symmetry operations 

172 nrot = len(self.rotations) 

173 retval.append(' %d symmetry operations (rot+trans)\n' % nrot) 

174 for i in range(nrot): 

175 retval.append(' ') 

176 for j in range(3): 

177 retval.append(' ') 

178 for k in range(3): 

179 retval.append(' %2d' % (self.rotations[i, j, k])) 

180 retval.append(' ') 

181 for j in range(3): 

182 retval.append(' %13.10f' % self.translations[i, j]) 

183 retval.append('\n') 

184 retval.append('\n') 

185 return ''.join(retval) 

186 

187 def __eq__(self, other): 

188 return self.no == other.no and self.setting == other.setting 

189 

190 def __ne__(self, other): 

191 return not self.__eq__(other) 

192 

193 def __lt__(self, other): 

194 return self.no < other.no or (self.no == other.no 

195 and self.setting < other.setting) 

196 

197 def __index__(self): 

198 return self.no 

199 

200 __int__ = __index__ 

201 

202 def get_symop(self): 

203 """Returns all symmetry operations (including inversions and 

204 subtranslations) as a sequence of (rotation, translation) 

205 tuples.""" 

206 symop = [] 

207 parities = [1] 

208 if self.centrosymmetric: 

209 parities.append(-1) 

210 for parity in parities: 

211 for subtrans in self.subtrans: 

212 for rot, trans in zip(self.rotations, self.translations): 

213 newtrans = np.mod(trans + subtrans, 1) 

214 symop.append((parity * rot, newtrans)) 

215 return symop 

216 

217 def get_op(self): 

218 """Returns all symmetry operations (including inversions and 

219 subtranslations), but unlike get_symop(), they are returned as 

220 two ndarrays.""" 

221 if self.centrosymmetric: 

222 rot = np.tile(np.vstack((self.rotations, -self.rotations)), 

223 (self.nsubtrans, 1, 1)) 

224 trans = np.tile(np.vstack((self.translations, -self.translations)), 

225 (self.nsubtrans, 1)) 

226 trans += np.repeat(self.subtrans, 2 * len(self.rotations), axis=0) 

227 trans = np.mod(trans, 1) 

228 else: 

229 rot = np.tile(self.rotations, (self.nsubtrans, 1, 1)) 

230 trans = np.tile(self.translations, (self.nsubtrans, 1)) 

231 trans += np.repeat(self.subtrans, len(self.rotations), axis=0) 

232 trans = np.mod(trans, 1) 

233 return rot, trans 

234 

235 def get_rotations(self): 

236 """Return all rotations, including inversions for 

237 centrosymmetric crystals.""" 

238 if self.centrosymmetric: 

239 return np.vstack((self.rotations, -self.rotations)) 

240 else: 

241 return self.rotations 

242 

243 def equivalent_reflections(self, hkl): 

244 """Return all equivalent reflections to the list of Miller indices 

245 in hkl. 

246 

247 Example: 

248 

249 >>> from ase.spacegroup import Spacegroup 

250 >>> sg = Spacegroup(225) # fcc 

251 >>> sg.equivalent_reflections([[0, 0, 2]]) 

252 array([[ 0, 0, -2], 

253 [ 0, -2, 0], 

254 [-2, 0, 0], 

255 [ 2, 0, 0], 

256 [ 0, 2, 0], 

257 [ 0, 0, 2]]) 

258 """ 

259 hkl = np.array(hkl, dtype='int', ndmin=2) 

260 rot = self.get_rotations() 

261 n, nrot = len(hkl), len(rot) 

262 R = rot.transpose(0, 2, 1).reshape((3 * nrot, 3)).T 

263 refl = np.dot(hkl, R).reshape((n * nrot, 3)) 

264 ind = np.lexsort(refl.T) 

265 refl = refl[ind] 

266 diff = np.diff(refl, axis=0) 

267 mask = np.any(diff, axis=1) 

268 return np.vstack((refl[:-1][mask], refl[-1, :])) 

269 

270 def equivalent_lattice_points(self, uvw): 

271 """Return all lattice points equivalent to any of the lattice points 

272 in `uvw` with respect to rotations only. 

273 

274 Only equivalent lattice points that conserves the distance to 

275 origo are included in the output (making this a kind of real 

276 space version of the equivalent_reflections() method). 

277 

278 Example: 

279 

280 >>> from ase.spacegroup import Spacegroup 

281 >>> sg = Spacegroup(225) # fcc 

282 >>> sg.equivalent_lattice_points([[0, 0, 2]]) 

283 array([[ 0, 0, -2], 

284 [ 0, -2, 0], 

285 [-2, 0, 0], 

286 [ 2, 0, 0], 

287 [ 0, 2, 0], 

288 [ 0, 0, 2]]) 

289 

290 """ 

291 uvw = np.array(uvw, ndmin=2) 

292 rot = self.get_rotations() 

293 n, nrot = len(uvw), len(rot) 

294 directions = np.dot(uvw, rot).reshape((n * nrot, 3)) 

295 ind = np.lexsort(directions.T) 

296 directions = directions[ind] 

297 diff = np.diff(directions, axis=0) 

298 mask = np.any(diff, axis=1) 

299 return np.vstack((directions[:-1][mask], directions[-1:])) 

300 

301 def symmetry_normalised_reflections(self, hkl): 

302 """Returns an array of same size as *hkl*, containing the 

303 corresponding symmetry-equivalent reflections of lowest 

304 indices. 

305 

306 Example: 

307 

308 >>> from ase.spacegroup import Spacegroup 

309 >>> sg = Spacegroup(225) # fcc 

310 >>> sg.symmetry_normalised_reflections([[2, 0, 0], [0, 2, 0]]) 

311 array([[ 0, 0, -2], 

312 [ 0, 0, -2]]) 

313 """ 

314 hkl = np.array(hkl, dtype=int, ndmin=2) 

315 normalised = np.empty(hkl.shape, int) 

316 R = self.get_rotations().transpose(0, 2, 1) 

317 for i, g in enumerate(hkl): 

318 gsym = np.dot(R, g) 

319 j = np.lexsort(gsym.T)[0] 

320 normalised[i, :] = gsym[j] 

321 return normalised 

322 

323 def unique_reflections(self, hkl): 

324 """Returns a subset *hkl* containing only the symmetry-unique 

325 reflections. 

326 

327 Example: 

328 

329 >>> from ase.spacegroup import Spacegroup 

330 >>> sg = Spacegroup(225) # fcc 

331 >>> sg.unique_reflections([[ 2, 0, 0], 

332 ... [ 0, -2, 0], 

333 ... [ 2, 2, 0], 

334 ... [ 0, -2, -2]]) 

335 array([[2, 0, 0], 

336 [2, 2, 0]]) 

337 """ 

338 hkl = np.array(hkl, dtype=int, ndmin=2) 

339 hklnorm = self.symmetry_normalised_reflections(hkl) 

340 perm = np.lexsort(hklnorm.T) 

341 iperm = perm.argsort() 

342 xmask = np.abs(np.diff(hklnorm[perm], axis=0)).any(axis=1) 

343 mask = np.concatenate(([True], xmask)) 

344 imask = mask[iperm] 

345 return hkl[imask] 

346 

347 def equivalent_sites(self, 

348 scaled_positions, 

349 onduplicates='error', 

350 symprec=1e-3, 

351 occupancies=None): 

352 """Returns the scaled positions and all their equivalent sites. 

353 

354 Parameters: 

355 

356 scaled_positions: list | array 

357 List of non-equivalent sites given in unit cell coordinates. 

358 

359 occupancies: list | array, optional (default=None) 

360 List of occupancies corresponding to the respective sites. 

361 

362 onduplicates : 'keep' | 'replace' | 'warn' | 'error' 

363 Action if `scaled_positions` contain symmetry-equivalent 

364 positions of full occupancy: 

365 

366 'keep' 

367 ignore additional symmetry-equivalent positions 

368 'replace' 

369 replace 

370 'warn' 

371 like 'keep', but issue an UserWarning 

372 'error' 

373 raises a SpacegroupValueError 

374 

375 symprec: float 

376 Minimum "distance" betweed two sites in scaled coordinates 

377 before they are counted as the same site. 

378 

379 Returns: 

380 

381 sites: array 

382 A NumPy array of equivalent sites. 

383 kinds: list 

384 A list of integer indices specifying which input site is 

385 equivalent to the corresponding returned site. 

386 

387 Example: 

388 

389 >>> from ase.spacegroup import Spacegroup 

390 >>> sg = Spacegroup(225) # fcc 

391 >>> sites, kinds = sg.equivalent_sites([[0, 0, 0], [0.5, 0.0, 0.0]]) 

392 >>> sites 

393 array([[ 0. , 0. , 0. ], 

394 [ 0. , 0.5, 0.5], 

395 [ 0.5, 0. , 0.5], 

396 [ 0.5, 0.5, 0. ], 

397 [ 0.5, 0. , 0. ], 

398 [ 0. , 0.5, 0. ], 

399 [ 0. , 0. , 0.5], 

400 [ 0.5, 0.5, 0.5]]) 

401 >>> kinds 

402 [0, 0, 0, 0, 1, 1, 1, 1] 

403 """ 

404 kinds = [] 

405 sites = [] 

406 

407 scaled = np.array(scaled_positions, ndmin=2) 

408 

409 for kind, pos in enumerate(scaled): 

410 for rot, trans in self.get_symop(): 

411 site = np.mod(np.dot(rot, pos) + trans, 1.) 

412 if not sites: 

413 sites.append(site) 

414 kinds.append(kind) 

415 continue 

416 t = site - sites 

417 mask = np.all( 

418 (abs(t) < symprec) | (abs(abs(t) - 1.0) < symprec), axis=1) 

419 if np.any(mask): 

420 inds = np.argwhere(mask).flatten() 

421 for ind in inds: 

422 # then we would just add the same thing again -> skip 

423 if kinds[ind] == kind: 

424 pass 

425 elif onduplicates == 'keep': 

426 pass 

427 elif onduplicates == 'replace': 

428 kinds[ind] = kind 

429 elif onduplicates == 'warn': 

430 warnings.warn('scaled_positions %d and %d ' 

431 'are equivalent' % 

432 (kinds[ind], kind)) 

433 elif onduplicates == 'error': 

434 raise SpacegroupValueError( 

435 'scaled_positions %d and %d are equivalent' % 

436 (kinds[ind], kind)) 

437 else: 

438 raise SpacegroupValueError( 

439 'Argument "onduplicates" must be one of: ' 

440 '"keep", "replace", "warn" or "error".') 

441 else: 

442 sites.append(site) 

443 kinds.append(kind) 

444 

445 return np.array(sites), kinds 

446 

447 def symmetry_normalised_sites(self, 

448 scaled_positions, 

449 map_to_unitcell=True): 

450 """Returns an array of same size as *scaled_positions*, 

451 containing the corresponding symmetry-equivalent sites of 

452 lowest indices. 

453 

454 If *map_to_unitcell* is true, the returned positions are all 

455 mapped into the unit cell, i.e. lattice translations are 

456 included as symmetry operator. 

457 

458 Example: 

459 

460 >>> from ase.spacegroup import Spacegroup 

461 >>> sg = Spacegroup(225) # fcc 

462 >>> sg.symmetry_normalised_sites([[0.0, 0.5, 0.5], [1.0, 1.0, 0.0]]) 

463 array([[ 0., 0., 0.], 

464 [ 0., 0., 0.]]) 

465 """ 

466 scaled = np.array(scaled_positions, ndmin=2) 

467 normalised = np.empty(scaled.shape, float) 

468 rot, trans = self.get_op() 

469 for i, pos in enumerate(scaled): 

470 sympos = np.dot(rot, pos) + trans 

471 if map_to_unitcell: 

472 # Must be done twice, see the scaled_positions.py test 

473 sympos %= 1.0 

474 sympos %= 1.0 

475 j = np.lexsort(sympos.T)[0] 

476 normalised[i, :] = sympos[j] 

477 return normalised 

478 

479 def unique_sites(self, 

480 scaled_positions, 

481 symprec=1e-3, 

482 output_mask=False, 

483 map_to_unitcell=True): 

484 """Returns a subset of *scaled_positions* containing only the 

485 symmetry-unique positions. If *output_mask* is True, a boolean 

486 array masking the subset is also returned. 

487 

488 If *map_to_unitcell* is true, all sites are first mapped into 

489 the unit cell making e.g. [0, 0, 0] and [1, 0, 0] equivalent. 

490 

491 Example: 

492 

493 >>> from ase.spacegroup import Spacegroup 

494 >>> sg = Spacegroup(225) # fcc 

495 >>> sg.unique_sites([[0.0, 0.0, 0.0], 

496 ... [0.5, 0.5, 0.0], 

497 ... [1.0, 0.0, 0.0], 

498 ... [0.5, 0.0, 0.0]]) 

499 array([[ 0. , 0. , 0. ], 

500 [ 0.5, 0. , 0. ]]) 

501 """ 

502 scaled = np.array(scaled_positions, ndmin=2) 

503 symnorm = self.symmetry_normalised_sites(scaled, map_to_unitcell) 

504 perm = np.lexsort(symnorm.T) 

505 iperm = perm.argsort() 

506 xmask = np.abs(np.diff(symnorm[perm], axis=0)).max(axis=1) > symprec 

507 mask = np.concatenate(([True], xmask)) 

508 imask = mask[iperm] 

509 if output_mask: 

510 return scaled[imask], imask 

511 else: 

512 return scaled[imask] 

513 

514 def tag_sites(self, scaled_positions, symprec=1e-3): 

515 """Returns an integer array of the same length as *scaled_positions*, 

516 tagging all equivalent atoms with the same index. 

517 

518 Example: 

519 

520 >>> from ase.spacegroup import Spacegroup 

521 >>> sg = Spacegroup(225) # fcc 

522 >>> sg.tag_sites([[0.0, 0.0, 0.0], 

523 ... [0.5, 0.5, 0.0], 

524 ... [1.0, 0.0, 0.0], 

525 ... [0.5, 0.0, 0.0]]) 

526 array([0, 0, 0, 1]) 

527 """ 

528 scaled = np.array(scaled_positions, ndmin=2) 

529 scaled %= 1.0 

530 scaled %= 1.0 

531 tags = -np.ones((len(scaled), ), dtype=int) 

532 mask = np.ones((len(scaled), ), dtype=bool) 

533 rot, trans = self.get_op() 

534 i = 0 

535 while mask.any(): 

536 pos = scaled[mask][0] 

537 sympos = np.dot(rot, pos) + trans 

538 # Must be done twice, see the scaled_positions.py test 

539 sympos %= 1.0 

540 sympos %= 1.0 

541 m = ~np.all(np.any(np.abs(scaled[np.newaxis, :, :] - 

542 sympos[:, np.newaxis, :]) > symprec, 

543 axis=2), 

544 axis=0) 

545 assert not np.any((~mask) & m) 

546 tags[m] = i 

547 mask &= ~m 

548 i += 1 

549 return tags 

550 

551 

552def get_datafile(): 

553 """Return default path to datafile.""" 

554 return os.path.join(os.path.dirname(__file__), 'spacegroup.dat') 

555 

556 

557def format_symbol(symbol): 

558 """Returns well formatted Hermann-Mauguin symbol as extected by 

559 the database, by correcting the case and adding missing or 

560 removing dublicated spaces.""" 

561 fixed = [] 

562 s = symbol.strip() 

563 s = s[0].upper() + s[1:].lower() 

564 for c in s: 

565 if c.isalpha(): 

566 if len(fixed) and fixed[-1] == '/': 

567 fixed.append(c) 

568 else: 

569 fixed.append(' ' + c + ' ') 

570 elif c.isspace(): 

571 fixed.append(' ') 

572 elif c.isdigit(): 

573 fixed.append(c) 

574 elif c == '-': 

575 fixed.append(' ' + c) 

576 elif c == '/': 

577 fixed.append(c) 

578 s = ''.join(fixed).strip() 

579 return ' '.join(s.split()) 

580 

581 

582# Functions for parsing the database. They are moved outside the 

583# Spacegroup class in order to make it easier to later implement 

584# caching to avoid reading the database each time a new Spacegroup 

585# instance is created. 

586 

587 

588def _skip_to_blank(f, spacegroup, setting): 

589 """Read lines from f until a blank line is encountered.""" 

590 while True: 

591 line = f.readline() 

592 if not line: 

593 raise SpacegroupNotFoundError( 

594 f'invalid spacegroup `{spacegroup}`, setting `{setting}` not ' 

595 'found in data base') 

596 if not line.strip(): 

597 break 

598 

599 

600def _skip_to_nonblank(f, spacegroup, setting): 

601 """Read lines from f until a nonblank line not starting with a 

602 hash (#) is encountered and returns this and the next line.""" 

603 while True: 

604 line1 = f.readline() 

605 if not line1: 

606 raise SpacegroupNotFoundError( 

607 'invalid spacegroup %s, setting %i not found in data base' % 

608 (spacegroup, setting)) 

609 line1.strip() 

610 if line1 and not line1.startswith('#'): 

611 line2 = f.readline() 

612 break 

613 return line1, line2 

614 

615 

616def _read_datafile_entry(spg, no, symbol, setting, f): 

617 """Read space group data from f to spg.""" 

618 

619 floats = {'0.0': 0.0, '1.0': 1.0, '0': 0.0, '1': 1.0, '-1': -1.0} 

620 for n, d in [(1, 2), (1, 3), (2, 3), (1, 4), (3, 4), (1, 6), (5, 6)]: 

621 floats[f'{n}/{d}'] = n / d 

622 floats[f'-{n}/{d}'] = -n / d 

623 

624 spg._no = no 

625 spg._symbol = symbol.strip() 

626 spg._setting = setting 

627 spg._centrosymmetric = bool(int(f.readline().split()[1])) 

628 # primitive vectors 

629 f.readline() 

630 spg._scaled_primitive_cell = np.array( 

631 [[float(floats.get(s, s)) for s in f.readline().split()] 

632 for i in range(3)], 

633 dtype=float) 

634 # primitive reciprocal vectors 

635 f.readline() 

636 spg._reciprocal_cell = np.array([[int(i) for i in f.readline().split()] 

637 for i in range(3)], 

638 dtype=int) 

639 # subtranslations 

640 spg._nsubtrans = int(f.readline().split()[0]) 

641 spg._subtrans = np.array( 

642 [[float(floats.get(t, t)) for t in f.readline().split()] 

643 for i in range(spg._nsubtrans)], 

644 dtype=float) 

645 # symmetry operations 

646 nsym = int(f.readline().split()[0]) 

647 symop = np.array([[float(floats.get(s, s)) for s in f.readline().split()] 

648 for i in range(nsym)], 

649 dtype=float) 

650 spg._nsymop = nsym 

651 spg._rotations = np.array(symop[:, :9].reshape((nsym, 3, 3)), dtype=int) 

652 spg._translations = symop[:, 9:] 

653 

654 

655def _read_datafile(spg, spacegroup, setting, f): 

656 if isinstance(spacegroup, int): 

657 pass 

658 elif isinstance(spacegroup, str): 

659 spacegroup = ' '.join(spacegroup.strip().split()) 

660 compact_spacegroup = ''.join(spacegroup.split()) 

661 else: 

662 raise SpacegroupValueError('`spacegroup` must be of type int or str') 

663 while True: 

664 line1, line2 = _skip_to_nonblank(f, spacegroup, setting) 

665 _no, _symbol = line1.strip().split(None, 1) 

666 _symbol = format_symbol(_symbol) 

667 compact_symbol = ''.join(_symbol.split()) 

668 _setting = int(line2.strip().split()[1]) 

669 _no = int(_no) 

670 

671 condition = ( 

672 (isinstance(spacegroup, int) and _no == spacegroup 

673 and _setting == setting) 

674 or (isinstance(spacegroup, str) 

675 and compact_symbol == compact_spacegroup) and 

676 (setting is None or _setting == setting)) 

677 

678 if condition: 

679 _read_datafile_entry(spg, _no, _symbol, _setting, f) 

680 break 

681 else: 

682 _skip_to_blank(f, spacegroup, setting) 

683 

684 

685def parse_sitesym_element(element): 

686 """Parses one element from a single site symmetry in the form used 

687 by the International Tables. 

688 

689 Examples: 

690 

691 >>> parse_sitesym_element("x") 

692 ([(0, 1)], 0.0) 

693 >>> parse_sitesym_element("-1/2-y") 

694 ([(1, -1)], -0.5) 

695 >>> parse_sitesym_element("z+0.25") 

696 ([(2, 1)], 0.25) 

697 >>> parse_sitesym_element("x-z+0.5") 

698 ([(0, 1), (2, -1)], 0.5) 

699 

700 

701 

702 Parameters 

703 ---------- 

704 

705 element: str 

706 Site symmetry like "x" or "-y+1/4" or "0.5+z". 

707 

708 

709 Returns 

710 ------- 

711 

712 list[tuple[int, int]] 

713 Rotation information in the form '(index, sign)' where index is 

714 0 for "x", 1 for "y" and 2 for "z" and sign is '1' for a positive 

715 entry and '-1' for a negative entry. E.g. "x" is '(0, 1)' and 

716 "-z" is (2, -1). 

717 

718 float 

719 Translation information in fractional space. E.g. "-1/4" is 

720 '-0.25' and "1/2" is '0.5' and "0.75" is '0.75'. 

721 

722 

723 """ 

724 element = element.lower() 

725 is_positive = True 

726 is_frac = False 

727 sng_trans = None 

728 fst_trans = [] 

729 snd_trans = [] 

730 rot = [] 

731 

732 for char in element: 

733 if char == "+": 

734 is_positive = True 

735 elif char == "-": 

736 is_positive = False 

737 elif char == "/": 

738 is_frac = True 

739 elif char in "xyz": 

740 rot.append((ord(char) - ord("x"), 1 if is_positive else -1)) 

741 elif char.isdigit() or char == ".": 

742 if sng_trans is None: 

743 sng_trans = 1.0 if is_positive else -1.0 

744 if is_frac: 

745 snd_trans.append(char) 

746 else: 

747 fst_trans.append(char) 

748 

749 trans = 0.0 if not fst_trans else (sng_trans * float("".join(fst_trans))) 

750 if is_frac: 

751 trans /= float("".join(snd_trans)) 

752 

753 return rot, trans 

754 

755 

756def parse_sitesym_single(sym, out_rot, out_trans, sep=",", 

757 force_positive_translation=False): 

758 """Parses a single site symmetry in the form used by International 

759 Tables and overwrites 'out_rot' and 'out_trans' with data. 

760 

761 Parameters 

762 ---------- 

763 

764 sym: str 

765 Site symmetry in the form used by International Tables 

766 (e.g. "x,y,z", "y-1/2,x,-z"). 

767 

768 out_rot: np.array 

769 A 3x3-integer array representing rotations (changes are made inplace). 

770 

771 out_rot: np.array 

772 A 3-float array representing translations (changes are made inplace). 

773 

774 sep: str 

775 String separator ("," in "x,y,z"). 

776 

777 force_positive_translation: bool 

778 Forces fractional translations to be between 0 and 1 (otherwise 

779 negative values might be accepted). Defaults to 'False'. 

780 

781 

782 Returns 

783 ------- 

784 

785 Nothing is returned: 'out_rot' and 'out_trans' are changed inplace. 

786 

787 

788 """ 

789 out_rot[:] = 0.0 

790 out_trans[:] = 0.0 

791 

792 for i, element in enumerate(sym.split(sep)): 

793 e_rot_list, e_trans = parse_sitesym_element(element) 

794 for rot_idx, rot_sgn in e_rot_list: 

795 out_rot[i][rot_idx] = rot_sgn 

796 out_trans[i] = \ 

797 (e_trans % 1.0) if force_positive_translation else e_trans 

798 

799 

800def parse_sitesym(symlist, sep=',', force_positive_translation=False): 

801 """Parses a sequence of site symmetries in the form used by 

802 International Tables and returns corresponding rotation and 

803 translation arrays. 

804 

805 Example: 

806 

807 >>> symlist = [ 

808 ... 'x,y,z', 

809 ... '-y+1/2,x+1/2,z', 

810 ... '-y,-x,-z', 

811 ... 'x-1/4, y-1/4, -z' 

812 ... ] 

813 >>> rot, trans = parse_sitesym(symlist) 

814 >>> rot 

815 array([[[ 1, 0, 0], 

816 [ 0, 1, 0], 

817 [ 0, 0, 1]], 

818 <BLANKLINE> 

819 [[ 0, -1, 0], 

820 [ 1, 0, 0], 

821 [ 0, 0, 1]], 

822 <BLANKLINE> 

823 [[ 0, -1, 0], 

824 [-1, 0, 0], 

825 [ 0, 0, -1]], 

826 <BLANKLINE> 

827 [[ 1, 0, 0], 

828 [ 0, 1, 0], 

829 [ 0, 0, -1]]]) 

830 >>> trans 

831 array([[ 0. , 0. , 0. ], 

832 [ 0.5 , 0.5 , 0. ], 

833 [ 0. , 0. , 0. ], 

834 [-0.25, -0.25, 0. ]]) 

835 """ 

836 

837 nsym = len(symlist) 

838 rot = np.zeros((nsym, 3, 3), dtype='int') 

839 trans = np.zeros((nsym, 3)) 

840 

841 for i, sym in enumerate(symlist): 

842 parse_sitesym_single( 

843 sym, rot[i], trans[i], sep=sep,