Coverage for /builds/debichem-team/python-ase/ase/spacegroup/spacegroup.py: 91.98%

424 statements  

« prev     ^ index     » next       coverage.py v7.5.3, created at 2025-03-06 04:00 +0000

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 lru_cache, total_ordering 

11from types import SimpleNamespace 

12from typing import Union 

13 

14import numpy as np 

15 

16from ase.utils import deprecated 

17 

18__all__ = ['Spacegroup'] 

19 

20 

21class SpacegroupError(Exception): 

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

23 

24 

25class SpacegroupNotFoundError(SpacegroupError): 

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

27 

28 

29class SpacegroupValueError(SpacegroupError): 

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

31 

32 

33# Type alias 

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

35 

36 

37@total_ordering 

38class Spacegroup: 

39 """A space group class. 

40 

41 The instances of Spacegroup describes the symmetry operations for 

42 the given space group. 

43 

44 Example: 

45 

46 >>> from ase.spacegroup import Spacegroup 

47 >>> 

48 >>> sg = Spacegroup(225) 

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

50 Space group 225 F m -3 m 

51 >>> sg.scaled_primitive_cell 

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

53 [ 0.5, 0. , 0.5], 

54 [ 0.5, 0.5, 0. ]]) 

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

56 >>> sites 

57 array([[ 0. , 0. , 0. ], 

58 [ 0. , 0.5, 0.5], 

59 [ 0.5, 0. , 0.5], 

60 [ 0.5, 0.5, 0. ]]) 

61 """ 

62 @property 

63 def no(self): 

64 """Space group number in International Tables of Crystallography.""" 

65 return self._no 

66 

67 @property 

68 def symbol(self): 

69 """Hermann-Mauguin (or international) symbol for the space group.""" 

70 return self._symbol 

71 

72 @property 

73 def setting(self): 

74 """Space group setting. Either one or two.""" 

75 return self._setting 

76 

77 @property 

78 def lattice(self): 

79 """Lattice type. 

80 

81 P primitive 

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

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

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

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

86 """ 

87 return self._symbol[0] 

88 

89 @property 

90 def centrosymmetric(self): 

91 """Whether a center of symmetry exists.""" 

92 return self._centrosymmetric 

93 

94 @property 

95 def scaled_primitive_cell(self): 

96 """Primitive cell in scaled coordinates. 

97 

98 Matrix with the primitive vectors along the rows. 

99 """ 

100 return self._scaled_primitive_cell 

101 

102 @property 

103 def reciprocal_cell(self): 

104 """ 

105 

106 Tree Miller indices that span all kinematically non-forbidden 

107 reflections as a matrix with the Miller indices along the rows. 

108 """ 

109 return self._reciprocal_cell 

110 

111 @property 

112 def nsubtrans(self): 

113 """Number of cell-subtranslation vectors.""" 

114 return len(self._subtrans) 

115 

116 @property 

117 def nsymop(self): 

118 """Total number of symmetry operations.""" 

119 scale = 2 if self.centrosymmetric else 1 

120 return scale * len(self._rotations) * len(self._subtrans) 

121 

122 @property 

123 def subtrans(self): 

124 """Translations vectors belonging to cell-sub-translations.""" 

125 return self._subtrans 

126 

127 @property 

128 def rotations(self): 

129 """Symmetry rotation matrices. 

130 

131 The invertions are not included for centrosymmetrical crystals. 

132 """ 

133 return self._rotations 

134 

135 @property 

136 def translations(self): 

137 """Symmetry translations. 

138 

139 The invertions are not included for centrosymmetrical crystals. 

140 """ 

141 return self._translations 

142 

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

144 """Returns a new Spacegroup instance. 

145 

146 Parameters: 

147 

148 spacegroup : int | string | Spacegroup instance 

149 The space group number in International Tables of 

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

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

152 setting : 1 | 2 

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

154 determines Which of these should be used. 

155 datafile : None | string 

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

157 will be used. 

158 """ 

159 if isinstance(spacegroup, Spacegroup): 

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

161 setattr(self, k, v) 

162 return 

163 if not datafile: 

164 datafile = get_datafile() 

165 namespace = _read_datafile(spacegroup, setting, datafile) 

166 self._no = namespace._no 

167 self._symbol = namespace._symbol 

168 self._setting = namespace._setting 

169 self._centrosymmetric = namespace._centrosymmetric 

170 self._scaled_primitive_cell = namespace._scaled_primitive_cell 

171 self._reciprocal_cell = namespace._reciprocal_cell 

172 self._subtrans = namespace._subtrans 

173 self._rotations = namespace._rotations 

174 self._translations = namespace._translations 

175 

176 def __repr__(self): 

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

178 

179 def todict(self): 

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

181 

182 def __str__(self): 

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

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

185 retval = [] 

186 # no, symbol 

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

188 # setting 

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

190 # centrosymmetric 

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

192 # primitive vectors 

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

194 for i in range(3): 

195 retval.append(' ') 

196 for j in range(3): 

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

198 retval.append('\n') 

199 # primitive reciprocal vectors 

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

201 for i in range(3): 

202 retval.append(' ') 

203 for j in range(3): 

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

205 retval.append('\n') 

206 # sublattice 

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

208 for i in range(self.nsubtrans): 

209 retval.append(' ') 

210 for j in range(3): 

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

212 retval.append('\n') 

213 # symmetry operations 

214 nrot = len(self.rotations) 

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

216 for i in range(nrot): 

217 retval.append(' ') 

218 for j in range(3): 

219 retval.append(' ') 

220 for k in range(3): 

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

222 retval.append(' ') 

223 for j in range(3): 

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

225 retval.append('\n') 

226 retval.append('\n') 

227 return ''.join(retval) 

228 

229 def __eq__(self, other): 

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

231 

232 def __ne__(self, other): 

233 return not self.__eq__(other) 

234 

235 def __lt__(self, other): 

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

237 and self.setting < other.setting) 

238 

239 def __index__(self): 

240 return self.no 

241 

242 __int__ = __index__ 

243 

244 def get_symop(self): 

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

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

247 tuples.""" 

248 symop = [] 

249 parities = [1] 

250 if self.centrosymmetric: 

251 parities.append(-1) 

252 for parity in parities: 

253 for subtrans in self.subtrans: 

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

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

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

257 return symop 

258 

259 def get_op(self): 

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

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

262 two ndarrays.""" 

263 if self.centrosymmetric: 

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

265 (self.nsubtrans, 1, 1)) 

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

267 (self.nsubtrans, 1)) 

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

269 trans = np.mod(trans, 1) 

270 else: 

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

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

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

274 trans = np.mod(trans, 1) 

275 return rot, trans 

276 

277 def get_rotations(self): 

278 """Return all rotations, including inversions for 

279 centrosymmetric crystals.""" 

280 if self.centrosymmetric: 

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

282 else: 

283 return self.rotations 

284 

285 def equivalent_reflections(self, hkl): 

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

287 in hkl. 

288 

289 Example: 

290 

291 >>> from ase.spacegroup import Spacegroup 

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

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

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

295 [ 0, -2, 0], 

296 [-2, 0, 0], 

297 [ 2, 0, 0], 

298 [ 0, 2, 0], 

299 [ 0, 0, 2]]) 

300 """ 

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

302 rot = self.get_rotations() 

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

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

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

306 ind = np.lexsort(refl.T) 

307 refl = refl[ind] 

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

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

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

311 

312 def equivalent_lattice_points(self, uvw): 

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

314 in `uvw` with respect to rotations only. 

315 

316 Only equivalent lattice points that conserves the distance to 

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

318 space version of the equivalent_reflections() method). 

319 

320 Example: 

321 

322 >>> from ase.spacegroup import Spacegroup 

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

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

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

326 [ 0, -2, 0], 

327 [-2, 0, 0], 

328 [ 2, 0, 0], 

329 [ 0, 2, 0], 

330 [ 0, 0, 2]]) 

331 

332 """ 

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

334 rot = self.get_rotations() 

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

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

337 ind = np.lexsort(directions.T) 

338 directions = directions[ind] 

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

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

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

342 

343 def symmetry_normalised_reflections(self, hkl): 

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

345 corresponding symmetry-equivalent reflections of lowest 

346 indices. 

347 

348 Example: 

349 

350 >>> from ase.spacegroup import Spacegroup 

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

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

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

354 [ 0, 0, -2]]) 

355 """ 

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

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

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

359 for i, g in enumerate(hkl): 

360 gsym = np.dot(R, g) 

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

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

363 return normalised 

364 

365 def unique_reflections(self, hkl): 

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

367 reflections. 

368 

369 Example: 

370 

371 >>> from ase.spacegroup import Spacegroup 

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

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

374 ... [ 0, -2, 0], 

375 ... [ 2, 2, 0], 

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

377 array([[2, 0, 0], 

378 [2, 2, 0]]) 

379 """ 

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

381 hklnorm = self.symmetry_normalised_reflections(hkl) 

382 perm = np.lexsort(hklnorm.T) 

383 iperm = perm.argsort() 

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

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

386 imask = mask[iperm] 

387 return hkl[imask] 

388 

389 def equivalent_sites(self, 

390 scaled_positions, 

391 onduplicates='error', 

392 symprec=1e-3, 

393 occupancies=None): 

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

395 

396 Parameters: 

397 

398 scaled_positions: list | array 

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

400 

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

402 List of occupancies corresponding to the respective sites. 

403 

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

405 Action if `scaled_positions` contain symmetry-equivalent 

406 positions of full occupancy: 

407 

408 'keep' 

409 ignore additional symmetry-equivalent positions 

410 'replace' 

411 replace 

412 'warn' 

413 like 'keep', but issue an UserWarning 

414 'error' 

415 raises a SpacegroupValueError 

416 

417 symprec: float 

418 Minimum "distance" betweed two sites in scaled coordinates 

419 before they are counted as the same site. 

420 

421 Returns: 

422 

423 sites: array 

424 A NumPy array of equivalent sites. 

425 kinds: list 

426 A list of integer indices specifying which input site is 

427 equivalent to the corresponding returned site. 

428 

429 Example: 

430 

431 >>> from ase.spacegroup import Spacegroup 

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

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

434 >>> sites 

435 array([[ 0. , 0. , 0. ], 

436 [ 0. , 0.5, 0.5], 

437 [ 0.5, 0. , 0.5], 

438 [ 0.5, 0.5, 0. ], 

439 [ 0.5, 0. , 0. ], 

440 [ 0. , 0.5, 0. ], 

441 [ 0. , 0. , 0.5], 

442 [ 0.5, 0.5, 0.5]]) 

443 >>> kinds 

444 [0, 0, 0, 0, 1, 1, 1, 1] 

445 """ 

446 if onduplicates not in ('keep', 'replace', 'warn', 'error'): 

447 raise SpacegroupValueError( 

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

449 '"keep", "replace", "warn" or "error".' 

450 ) 

451 

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

453 rotations, translations = zip(*self.get_symop()) 

454 rotations = np.array(rotations) 

455 translations = np.array(translations) 

456 

457 def find_orbit(point: np.ndarray) -> np.ndarray: 

458 """Find crystallographic orbit of the given point.""" 

459 candidates = ((rotations @ point) + translations % 1.0) % 1.0 

460 orbit = [candidates[0]] 

461 for member in candidates[1:]: 

462 diff = member - orbit 

463 diff -= np.rint(diff) 

464 if not np.any(np.all(np.abs(diff) < symprec, axis=1)): 

465 orbit.append(member) 

466 return np.array(orbit) 

467 

468 orbits = [] 

469 for kind, pos in enumerate(scaled): 

470 for i, (kind0, positions0) in enumerate(orbits): 

471 diff = pos - positions0 

472 diff -= np.rint(diff) 

473 if np.any(np.all(np.abs(diff) < symprec, axis=1)): 

474 if onduplicates == 'keep': 

475 pass 

476 elif onduplicates == 'replace': 

477 orbits[i] = (kind, positions0) 

478 elif onduplicates == 'warn': 

479 warnings.warn( 

480 'scaled_positions %d and %d are equivalent' % 

481 (kind0, kind)) 

482 elif onduplicates == 'error': 

483 raise SpacegroupValueError( 

484 'scaled_positions %d and %d are equivalent' % 

485 (kind0, kind)) 

486 break 

487 else: 

488 orbits.append((kind, find_orbit(pos))) 

489 

490 kinds = [] 

491 sites = [] 

492 for kind, orbit in orbits: 

493 kinds.extend(len(orbit) * [kind]) 

494 sites.append(orbit) 

495 

496 return np.concatenate(sites, axis=0), kinds 

497 

498 def symmetry_normalised_sites(self, 

499 scaled_positions, 

500 map_to_unitcell=True): 

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

502 containing the corresponding symmetry-equivalent sites of 

503 lowest indices. 

504 

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

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

507 included as symmetry operator. 

508 

509 Example: 

510 

511 >>> from ase.spacegroup import Spacegroup 

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

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

514 array([[ 0., 0., 0.], 

515 [ 0., 0., 0.]]) 

516 """ 

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

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

519 rot, trans = self.get_op() 

520 for i, pos in enumerate(scaled): 

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

522 if map_to_unitcell: 

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

524 sympos %= 1.0 

525 sympos %= 1.0 

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

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

528 return normalised 

529 

530 def unique_sites(self, 

531 scaled_positions, 

532 symprec=1e-3, 

533 output_mask=False, 

534 map_to_unitcell=True): 

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

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

537 array masking the subset is also returned. 

538 

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

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

541 

542 Example: 

543 

544 >>> from ase.spacegroup import Spacegroup 

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

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

547 ... [0.5, 0.5, 0.0], 

548 ... [1.0, 0.0, 0.0], 

549 ... [0.5, 0.0, 0.0]]) 

550 array([[ 0. , 0. , 0. ], 

551 [ 0.5, 0. , 0. ]]) 

552 """ 

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

554 symnorm = self.symmetry_normalised_sites(scaled, map_to_unitcell) 

555 perm = np.lexsort(symnorm.T) 

556 iperm = perm.argsort() 

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

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

559 imask = mask[iperm] 

560 if output_mask: 

561 return scaled[imask], imask 

562 else: 

563 return scaled[imask] 

564 

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

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

567 tagging all equivalent atoms with the same index. 

568 

569 Example: 

570 

571 >>> from ase.spacegroup import Spacegroup 

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

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

574 ... [0.5, 0.5, 0.0], 

575 ... [1.0, 0.0, 0.0], 

576 ... [0.5, 0.0, 0.0]]) 

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

578 """ 

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

580 scaled %= 1.0 

581 scaled %= 1.0 

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

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

584 rot, trans = self.get_op() 

585 i = 0 

586 while mask.any(): 

587 pos = scaled[mask][0] 

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

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

590 sympos %= 1.0 

591 sympos %= 1.0 

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

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

594 axis=2), 

595 axis=0) 

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

597 tags[m] = i 

598 mask &= ~m 

599 i += 1 

600 return tags 

601 

602 

603def get_datafile(): 

604 """Return default path to datafile.""" 

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

606 

607 

608def format_symbol(symbol): 

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

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

611 removing dublicated spaces.""" 

612 fixed = [] 

613 s = symbol.strip() 

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

615 for c in s: 

616 if c.isalpha(): 

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

618 fixed.append(c) 

619 else: 

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

621 elif c.isspace(): 

622 fixed.append(' ') 

623 elif c.isdigit(): 

624 fixed.append(c) 

625 elif c == '-': 

626 fixed.append(' ' + c) 

627 elif c == '/': 

628 fixed.append(c) 

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

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

631 

632 

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

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

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

636# instance is created. 

637 

638 

639def _skip_to_blank(f, spacegroup, setting): 

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

641 while True: 

642 line = f.readline() 

643 if not line: 

644 raise SpacegroupNotFoundError( 

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

646 'found in data base') 

647 if not line.strip(): 

648 break 

649 

650 

651def _skip_to_nonblank(f, spacegroup, setting): 

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

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

654 while True: 

655 line1 = f.readline() 

656 if not line1: 

657 raise SpacegroupNotFoundError( 

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

659 (spacegroup, setting)) 

660 line1.strip() 

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

662 line2 = f.readline() 

663 break 

664 return line1, line2 

665 

666 

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

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

669 

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

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

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

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

674 

675 spg._no = no 

676 spg._symbol = symbol.strip() 

677 spg._setting = setting 

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

679 # primitive vectors 

680 f.readline() 

681 spg._scaled_primitive_cell = np.array( 

682 [ 

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

684 for _ in range(3) 

685 ], 

686 dtype=float, 

687 ) 

688 # primitive reciprocal vectors 

689 f.readline() 

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

691 for i in range(3)], 

692 dtype=int) 

693 # subtranslations 

694 nsubtrans = int(f.readline().split()[0]) 

695 spg._subtrans = np.array( 

696 [ 

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

698 for _ in range(nsubtrans) 

699 ], 

700 dtype=float, 

701 ) 

702 # symmetry operations 

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

704 symop = np.array( 

705 [ 

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

707 for _ in range(nsym) 

708 ], 

709 dtype=float, 

710 ) 

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

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

713 

714 

715@lru_cache 

716def _read_datafile(spacegroup, setting, datafile): 

717 with open(datafile, encoding='utf-8') as fd: 

718 return _read_f(spacegroup, setting, fd) 

719 

720 

721def _read_f(spacegroup, setting, f): 

722 if isinstance(spacegroup, int): 

723 pass 

724 elif isinstance(spacegroup, str): 

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

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

727 else: 

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

729 while True: 

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

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

732 _symbol = format_symbol(_symbol) 

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

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

735 _no = int(_no) 

736 

737 condition = ( 

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

739 and _setting == setting) 

740 or (isinstance(spacegroup, str) 

741 and compact_symbol == compact_spacegroup) and 

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

743 

744 if condition: 

745 namespace = SimpleNamespace() 

746 _read_datafile_entry(namespace, _no, _symbol, _setting, f) 

747 return namespace 

748 else: 

749 _skip_to_blank(f, spacegroup, setting) 

750 

751 

752def parse_sitesym_element(element): 

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

754 by the International Tables. 

755 

756 Examples: 

757 

758 >>> parse_sitesym_element("x") 

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

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

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

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

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

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

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

766 

767 

768 

769 Parameters 

770 ---------- 

771 

772 element: str 

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

774 

775 

776 Returns 

777 ------- 

778 

779 list[tuple[int, int]] 

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

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

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

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

784 

785 float 

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

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

788 

789 

790 """ 

791 element = element.lower() 

792 is_positive = True 

793 is_frac = False 

794 sng_trans = None 

795 fst_trans = [] 

796 snd_trans = [] 

797 rot = [] 

798 

799 for char in element: 

800 if char == "+": 

801 is_positive = True 

802 elif char == "-": 

803 is_positive = False 

804 elif char == "/": 

805 is_frac = True 

806 elif char in "xyz": 

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

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

809 if sng_trans is None: 

810 sng_trans = 1.0 if is_positive else -1.0 

811 if is_frac: 

812 snd_trans.append(char) 

813 else: 

814 fst_trans.append(char) 

815 

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

817 if is_frac: 

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

819 

820 return rot, trans 

821 

822 

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

824 force_positive_translation=False): 

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

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

827 

828 Parameters 

829 ---------- 

830 

831 sym: str 

832 Site symmetry in the form used by International Tables 

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

834 

835 out_rot: np.array 

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

837 

838 out_rot: np.array 

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

840 

841 sep: str 

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

843 

844 force_positive_translation: bool 

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

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

847 

848 

849 Returns 

850 ------- 

851 

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

853 

854 

855 """ 

856 out_rot[:] = 0.0 

857 out_trans[:] = 0.0 

858 

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

860 e_rot_list, e_trans = parse_sitesym_element(element) 

861 for rot_idx, rot_sgn in e_rot_list: 

862 out_rot[i][rot_idx] = rot_sgn 

863 out_trans[i] = \ 

864 (e_trans % 1.0) if force_positive_translation else e_trans 

865 

866 

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

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

869 International Tables and returns corresponding rotation and 

870 translation arrays. 

871 

872 Example: 

873 

874 >>> symlist = [ 

875 ... 'x,y,z', 

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

877 ... '-y,-x,-z', 

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

879 ... ] 

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

881 >>> rot 

882 array([[[ 1, 0, 0], 

883 [ 0, 1, 0], 

884 [ 0, 0, 1]], 

885 <BLANKLINE> 

886 [[ 0, -1, 0], 

887 [ 1, 0, 0], 

888 [ 0, 0, 1]], 

889 <BLANKLINE> 

890 [[ 0, -1, 0], 

891 [-1, 0, 0], 

892 [ 0, 0, -1]], 

893 <BLANKLINE> 

894 [[ 1, 0, 0], 

895 [ 0, 1, 0], 

896 [ 0, 0, -1]]]) 

897 >>> trans 

898 array([[ 0. , 0. , 0. ], 

899 [ 0.5 , 0.5 , 0. ], 

900 [ 0. , 0. , 0. ], 

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

902 """ 

903 

904 nsym = len(symlist) 

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

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

907 

908 for i, sym in enumerate(symlist): 

909 parse_sitesym_single( 

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

911 force_positive_translation=force_positive_translation) 

912 

913 return rot, trans 

914 

915 

916def spacegroup_from_data(no=None, 

917 symbol=None, 

918 setting=None, 

919 centrosymmetric=None, 

920 scaled_primitive_cell=None, 

921 reciprocal_cell=None, 

922 subtrans=None, 

923 sitesym=None, 

924 rotations=None, 

925 translations=None, 

926 datafile=None): 

927 """Manually create a new space group instance. This might be 

928 useful when reading crystal data with its own spacegroup 

929 definitions.""" 

930 if no is not None and setting is not None: 

931 spg = Spacegroup(no, setting, datafile) 

932 elif symbol is not None: 

933 spg = Spacegroup(symbol, None, datafile) 

934 else: 

935 raise SpacegroupValueError('either *no* and *setting* ' 

936 'or *symbol* must be given') 

937 if not isinstance(sitesym, list): 

938 raise TypeError('sitesym must be a list') 

939 

940 have_sym = False 

941 if centrosymmetric is not None: 

942 spg._centrosymmetric = bool(centrosymmetric) 

943 if scaled_primitive_cell is not None: 

944 spg._scaled_primitive_cell = np.array(scaled_primitive_cell) 

945 if reciprocal_cell is not None: 

946 spg._reciprocal_cell = np.array(reciprocal_cell) 

947 if subtrans is not None: 

948 spg._subtrans = np.atleast_2d(subtrans) 

949 if sitesym is not None: 

950 spg._rotations, spg._translations = parse_sitesym(sitesym) 

951 have_sym = True 

952 if rotations is not None: 

953 spg._rotations = np.atleast_3d(rotations) 

954 have_sym = True 

955 if translations is not None: 

956 spg._translations = np.atleast_2d(translations) 

957 have_sym = True 

958 if have_sym: 

959 if spg._rotations.shape[0] != spg._translations.shape[0]: 

960 raise SpacegroupValueError('inconsistent number of rotations and ' 

961 'translations') 

962 return spg 

963 

964 

965@deprecated( 

966 '`get_spacegroup` has been deprecated due to its misleading output. ' 

967 'The returned `Spacegroup` object has symmetry operations for a ' 

968 'standard setting regardress of the given `Atoms` object. ' 

969 'See https://gitlab.com/ase/ase/-/issues/1534 for details. ' 

970 'Please use `ase.spacegroup.symmetrize.check_symmetry` or `spglib` ' 

971 'directly to get the symmetry operations for the given `Atoms` object.' 

972) 

973def get_spacegroup(atoms, symprec=1e-5): 

974 """Determine the spacegroup to which belongs the Atoms object. 

975 

976 This requires spglib: https://atztogo.github.io/spglib/ . 

977 

978 .. warning:: 

979 The returned ``Spacegroup`` object has symmetry operations for a 

980 standard setting regardless of the given ``Atoms`` object. 

981 See https://gitlab.com/ase/ase/-/issues/1534 for details. 

982 

983 .. deprecated:: 3.24.0 

984 Please use ``ase.spacegroup.symmetrize.check_symmetry`` or ``spglib`` 

985 directly to get the symmetry operations for the given ``Atoms`` object. 

986 

987 Parameters: 

988 

989 atoms: Atoms object 

990 Types, positions and unit-cell. 

991 symprec: float 

992 Symmetry tolerance, i.e. distance tolerance in Cartesian 

993 coordinates to find crystal symmetry. 

994 

995 The Spacegroup object is returned. 

996 """ 

997 

998 # Example: 

999 # (We don't include the example in docstring to appease doctests 

1000 # when import fails) 

1001 # >>> from ase.build import bulk 

1002 # >>> atoms = bulk("Cu", "fcc", a=3.6, cubic=True) 

1003 # >>> sg = get_spacegroup(atoms) 

1004 # >>> sg 

1005 # Spacegroup(225, setting=1) 

1006 # >>> sg.no 

1007 # 225 

1008 

1009 import spglib 

1010 

1011 sg = spglib.get_spacegroup((atoms.get_cell(), atoms.get_scaled_positions(), 

1012 atoms.get_atomic_numbers()), 

1013 symprec=symprec) 

1014 if sg is None: 

1015 raise RuntimeError('Spacegroup not found') 

1016 sg_no = int(sg[sg.find('(') + 1:sg.find(')')]) 

1017 return Spacegroup(sg_no)