Coverage for ase / spacegroup / spacegroup.py: 91.98%

424 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-03-30 08:22 +0000

1# fmt: off 

2 

3# Copyright (C) 2010, Jesper Friis 

4# (see accompanying license files for details). 

5"""Definition of the Spacegroup class. 

6 

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

8""" 

9 

10import os 

11import warnings 

12from functools import lru_cache, total_ordering 

13from types import SimpleNamespace 

14from typing import Union 

15 

16import numpy as np 

17 

18from ase.utils import deprecated, spglib_new_errorhandling 

19 

20__all__ = ['Spacegroup'] 

21 

22 

23class SpacegroupError(Exception): 

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

25 

26 

27class SpacegroupNotFoundError(SpacegroupError): 

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

29 

30 

31class SpacegroupValueError(SpacegroupError): 

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

33 

34 

35# Type alias 

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

37 

38 

39@total_ordering 

40class Spacegroup: 

41 """A space group class. 

42 

43 The instances of Spacegroup describes the symmetry operations for 

44 the given space group. 

45 

46 Example: 

47 

48 >>> from ase.spacegroup import Spacegroup 

49 >>> 

50 >>> sg = Spacegroup(225) 

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

52 Space group 225 F m -3 m 

53 >>> sg.scaled_primitive_cell 

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

55 [ 0.5, 0. , 0.5], 

56 [ 0.5, 0.5, 0. ]]) 

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

58 >>> sites 

59 array([[ 0. , 0. , 0. ], 

60 [ 0. , 0.5, 0.5], 

61 [ 0.5, 0. , 0.5], 

62 [ 0.5, 0.5, 0. ]]) 

63 """ 

64 @property 

65 def no(self): 

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

67 return self._no 

68 

69 @property 

70 def symbol(self): 

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

72 return self._symbol 

73 

74 @property 

75 def setting(self): 

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

77 return self._setting 

78 

79 @property 

80 def lattice(self): 

81 """Lattice type. 

82 

83 P primitive 

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

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

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

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

88 """ 

89 return self._symbol[0] 

90 

91 @property 

92 def centrosymmetric(self): 

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

94 return self._centrosymmetric 

95 

96 @property 

97 def scaled_primitive_cell(self): 

98 """Primitive cell in scaled coordinates. 

99 

100 Matrix with the primitive vectors along the rows. 

101 """ 

102 return self._scaled_primitive_cell 

103 

104 @property 

105 def reciprocal_cell(self): 

106 """ 

107 

108 Tree Miller indices that span all kinematically non-forbidden 

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

110 """ 

111 return self._reciprocal_cell 

112 

113 @property 

114 def nsubtrans(self): 

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

116 return len(self._subtrans) 

117 

118 @property 

119 def nsymop(self): 

120 """Total number of symmetry operations.""" 

121 scale = 2 if self.centrosymmetric else 1 

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

123 

124 @property 

125 def subtrans(self): 

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

127 return self._subtrans 

128 

129 @property 

130 def rotations(self): 

131 """Symmetry rotation matrices. 

132 

133 The invertions are not included for centrosymmetrical crystals. 

134 """ 

135 return self._rotations 

136 

137 @property 

138 def translations(self): 

139 """Symmetry translations. 

140 

141 The invertions are not included for centrosymmetrical crystals. 

142 """ 

143 return self._translations 

144 

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

146 """Returns a new Spacegroup instance. 

147 

148 Parameters 

149 ---------- 

150 

151 spacegroup : int | string | Spacegroup instance 

152 The space group number in International Tables of 

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

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

155 setting : 1 | 2 

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

157 determines Which of these should be used. 

158 datafile : None | string 

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

160 will be used. 

161 """ 

162 if isinstance(spacegroup, Spacegroup): 

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

164 setattr(self, k, v) 

165 return 

166 if not datafile: 

167 datafile = get_datafile() 

168 namespace = _read_datafile(spacegroup, setting, datafile) 

169 self._no = namespace._no 

170 self._symbol = namespace._symbol 

171 self._setting = namespace._setting 

172 self._centrosymmetric = namespace._centrosymmetric 

173 self._scaled_primitive_cell = namespace._scaled_primitive_cell 

174 self._reciprocal_cell = namespace._reciprocal_cell 

175 self._subtrans = namespace._subtrans 

176 self._rotations = namespace._rotations 

177 self._translations = namespace._translations 

178 

179 def __repr__(self): 

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

181 

182 def todict(self): 

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

184 

185 def __str__(self): 

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

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

188 retval = [] 

189 # no, symbol 

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

191 # setting 

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

193 # centrosymmetric 

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

195 # primitive vectors 

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

197 for i in range(3): 

198 retval.append(' ') 

199 for j in range(3): 

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

201 retval.append('\n') 

202 # primitive reciprocal vectors 

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

204 for i in range(3): 

205 retval.append(' ') 

206 for j in range(3): 

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

208 retval.append('\n') 

209 # sublattice 

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

211 for i in range(self.nsubtrans): 

212 retval.append(' ') 

213 for j in range(3): 

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

215 retval.append('\n') 

216 # symmetry operations 

217 nrot = len(self.rotations) 

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

219 for i in range(nrot): 

220 retval.append(' ') 

221 for j in range(3): 

222 retval.append(' ') 

223 for k in range(3): 

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

225 retval.append(' ') 

226 for j in range(3): 

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

228 retval.append('\n') 

229 retval.append('\n') 

230 return ''.join(retval) 

231 

232 def __eq__(self, other): 

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

234 

235 def __ne__(self, other): 

236 return not self.__eq__(other) 

237 

238 def __lt__(self, other): 

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

240 and self.setting < other.setting) 

241 

242 def __index__(self): 

243 return self.no 

244 

245 __int__ = __index__ 

246 

247 def get_symop(self): 

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

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

250 tuples.""" 

251 symop = [] 

252 parities = [1] 

253 if self.centrosymmetric: 

254 parities.append(-1) 

255 for parity in parities: 

256 for subtrans in self.subtrans: 

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

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

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

260 return symop 

261 

262 def get_op(self): 

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

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

265 two ndarrays.""" 

266 if self.centrosymmetric: 

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

268 (self.nsubtrans, 1, 1)) 

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

270 (self.nsubtrans, 1)) 

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

272 trans = np.mod(trans, 1) 

273 else: 

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

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

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

277 trans = np.mod(trans, 1) 

278 return rot, trans 

279 

280 def get_rotations(self): 

281 """Return all rotations, including inversions for 

282 centrosymmetric crystals.""" 

283 if self.centrosymmetric: 

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

285 else: 

286 return self.rotations 

287 

288 def equivalent_reflections(self, hkl): 

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

290 in hkl. 

291 

292 Example: 

293 

294 >>> from ase.spacegroup import Spacegroup 

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

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

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

298 [ 0, -2, 0], 

299 [-2, 0, 0], 

300 [ 2, 0, 0], 

301 [ 0, 2, 0], 

302 [ 0, 0, 2]]) 

303 """ 

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

305 rot = self.get_rotations() 

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

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

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

309 ind = np.lexsort(refl.T) 

310 refl = refl[ind] 

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

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

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

314 

315 def equivalent_lattice_points(self, uvw): 

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

317 in `uvw` with respect to rotations only. 

318 

319 Only equivalent lattice points that conserves the distance to 

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

321 space version of the equivalent_reflections() method). 

322 

323 Example: 

324 

325 >>> from ase.spacegroup import Spacegroup 

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

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

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

329 [ 0, -2, 0], 

330 [-2, 0, 0], 

331 [ 2, 0, 0], 

332 [ 0, 2, 0], 

333 [ 0, 0, 2]]) 

334 

335 """ 

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

337 rot = self.get_rotations() 

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

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

340 ind = np.lexsort(directions.T) 

341 directions = directions[ind] 

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

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

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

345 

346 def symmetry_normalised_reflections(self, hkl): 

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

348 corresponding symmetry-equivalent reflections of lowest 

349 indices. 

350 

351 Example: 

352 

353 >>> from ase.spacegroup import Spacegroup 

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

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

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

357 [ 0, 0, -2]]) 

358 """ 

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

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

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

362 for i, g in enumerate(hkl): 

363 gsym = np.dot(R, g) 

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

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

366 return normalised 

367 

368 def unique_reflections(self, hkl): 

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

370 reflections. 

371 

372 Example: 

373 

374 >>> from ase.spacegroup import Spacegroup 

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

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

377 ... [ 0, -2, 0], 

378 ... [ 2, 2, 0], 

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

380 array([[2, 0, 0], 

381 [2, 2, 0]]) 

382 """ 

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

384 hklnorm = self.symmetry_normalised_reflections(hkl) 

385 perm = np.lexsort(hklnorm.T) 

386 iperm = perm.argsort() 

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

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

389 imask = mask[iperm] 

390 return hkl[imask] 

391 

392 def equivalent_sites(self, 

393 scaled_positions, 

394 onduplicates='error', 

395 symprec=1e-3, 

396 occupancies=None): 

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

398 

399 Parameters 

400 ---------- 

401 

402 scaled_positions: list | array 

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

404 

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

406 List of occupancies corresponding to the respective sites. 

407 

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

409 Action if `scaled_positions` contain symmetry-equivalent 

410 positions of full occupancy: 

411 

412 'keep' 

413 ignore additional symmetry-equivalent positions 

414 'replace' 

415 replace 

416 'warn' 

417 like 'keep', but issue an UserWarning 

418 'error' 

419 raises a SpacegroupValueError 

420 

421 symprec: float 

422 Minimum "distance" betweed two sites in scaled coordinates 

423 before they are counted as the same site. 

424 

425 Returns 

426 ------- 

427 

428 sites: array 

429 A NumPy array of equivalent sites. 

430 kinds: list 

431 A list of integer indices specifying which input site is 

432 equivalent to the corresponding returned site. 

433 

434 Example: 

435 

436 >>> from ase.spacegroup import Spacegroup 

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

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

439 >>> sites 

440 array([[ 0. , 0. , 0. ], 

441 [ 0. , 0.5, 0.5], 

442 [ 0.5, 0. , 0.5], 

443 [ 0.5, 0.5, 0. ], 

444 [ 0.5, 0. , 0. ], 

445 [ 0. , 0.5, 0. ], 

446 [ 0. , 0. , 0.5], 

447 [ 0.5, 0.5, 0.5]]) 

448 >>> kinds 

449 [0, 0, 0, 0, 1, 1, 1, 1] 

450 """ 

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

452 raise SpacegroupValueError( 

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

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

455 ) 

456 

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

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

459 rotations = np.array(rotations) 

460 translations = np.array(translations) 

461 

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

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

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

465 orbit = [candidates[0]] 

466 for member in candidates[1:]: 

467 diff = member - orbit 

468 diff -= np.rint(diff) 

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

470 orbit.append(member) 

471 return np.array(orbit) 

472 

473 orbits = [] 

474 for kind, pos in enumerate(scaled): 

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

476 diff = pos - positions0 

477 diff -= np.rint(diff) 

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

479 if onduplicates == 'keep': 

480 pass 

481 elif onduplicates == 'replace': 

482 orbits[i] = (kind, positions0) 

483 elif onduplicates == 'warn': 

484 warnings.warn( 

485 'scaled_positions %d and %d are equivalent' % 

486 (kind0, kind)) 

487 elif onduplicates == 'error': 

488 raise SpacegroupValueError( 

489 'scaled_positions %d and %d are equivalent' % 

490 (kind0, kind)) 

491 break 

492 else: 

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

494 

495 kinds = [] 

496 sites = [] 

497 for kind, orbit in orbits: 

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

499 sites.append(orbit) 

500 

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

502 

503 def symmetry_normalised_sites(self, 

504 scaled_positions, 

505 map_to_unitcell=True): 

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

507 containing the corresponding symmetry-equivalent sites of 

508 lowest indices. 

509 

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

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

512 included as symmetry operator. 

513 

514 Example: 

515 

516 >>> from ase.spacegroup import Spacegroup 

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

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

519 array([[ 0., 0., 0.], 

520 [ 0., 0., 0.]]) 

521 """ 

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

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

524 rot, trans = self.get_op() 

525 for i, pos in enumerate(scaled): 

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

527 if map_to_unitcell: 

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

529 sympos %= 1.0 

530 sympos %= 1.0 

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

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

533 return normalised 

534 

535 def unique_sites(self, 

536 scaled_positions, 

537 symprec=1e-3, 

538 output_mask=False, 

539 map_to_unitcell=True): 

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

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

542 array masking the subset is also returned. 

543 

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

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

546 

547 Example: 

548 

549 >>> from ase.spacegroup import Spacegroup 

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

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

552 ... [0.5, 0.5, 0.0], 

553 ... [1.0, 0.0, 0.0], 

554 ... [0.5, 0.0, 0.0]]) 

555 array([[ 0. , 0. , 0. ], 

556 [ 0.5, 0. , 0. ]]) 

557 """ 

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

559 symnorm = self.symmetry_normalised_sites(scaled, map_to_unitcell) 

560 perm = np.lexsort(symnorm.T) 

561 iperm = perm.argsort() 

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

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

564 imask = mask[iperm] 

565 if output_mask: 

566 return scaled[imask], imask 

567 else: 

568 return scaled[imask] 

569 

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

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

572 tagging all equivalent atoms with the same index. 

573 

574 Example: 

575 

576 >>> from ase.spacegroup import Spacegroup 

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

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

579 ... [0.5, 0.5, 0.0], 

580 ... [1.0, 0.0, 0.0], 

581 ... [0.5, 0.0, 0.0]]) 

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

583 """ 

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

585 scaled %= 1.0 

586 scaled %= 1.0 

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

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

589 rot, trans = self.get_op() 

590 i = 0 

591 while mask.any(): 

592 pos = scaled[mask][0] 

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

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

595 sympos %= 1.0 

596 sympos %= 1.0 

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

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

599 axis=2), 

600 axis=0) 

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

602 tags[m] = i 

603 mask &= ~m 

604 i += 1 

605 return tags 

606 

607 

608def get_datafile(): 

609 """Return default path to datafile.""" 

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

611 

612 

613def format_symbol(symbol): 

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

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

616 removing dublicated spaces.""" 

617 fixed = [] 

618 s = symbol.strip() 

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

620 for c in s: 

621 if c.isalpha(): 

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

623 fixed.append(c) 

624 else: 

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

626 elif c.isspace(): 

627 fixed.append(' ') 

628 elif c.isdigit(): 

629 fixed.append(c) 

630 elif c == '-': 

631 fixed.append(' ' + c) 

632 elif c == '/': 

633 fixed.append(c) 

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

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

636 

637 

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

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

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

641# instance is created. 

642 

643 

644def _skip_to_blank(f, spacegroup, setting): 

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

646 while True: 

647 line = f.readline() 

648 if not line: 

649 raise SpacegroupNotFoundError( 

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

651 'found in data base') 

652 if not line.strip(): 

653 break 

654 

655 

656def _skip_to_nonblank(f, spacegroup, setting): 

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

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

659 while True: 

660 line1 = f.readline() 

661 if not line1: 

662 raise SpacegroupNotFoundError( 

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

664 (spacegroup, setting)) 

665 line1.strip() 

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

667 line2 = f.readline() 

668 break 

669 return line1, line2 

670 

671 

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

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

674 

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

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

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

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

679 

680 spg._no = no 

681 spg._symbol = symbol.strip() 

682 spg._setting = setting 

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

684 # primitive vectors 

685 f.readline() 

686 spg._scaled_primitive_cell = np.array( 

687 [ 

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

689 for _ in range(3) 

690 ], 

691 dtype=float, 

692 ) 

693 # primitive reciprocal vectors 

694 f.readline() 

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

696 for i in range(3)], 

697 dtype=int) 

698 # subtranslations 

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

700 spg._subtrans = np.array( 

701 [ 

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

703 for _ in range(nsubtrans) 

704 ], 

705 dtype=float, 

706 ) 

707 # symmetry operations 

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

709 symop = np.array( 

710 [ 

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

712 for _ in range(nsym) 

713 ], 

714 dtype=float, 

715 ) 

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

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

718 

719 

720@lru_cache 

721def _read_datafile(spacegroup, setting, datafile): 

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

723 return _read_f(spacegroup, setting, fd) 

724 

725 

726def _read_f(spacegroup, setting, f): 

727 if isinstance(spacegroup, int): 

728 pass 

729 elif isinstance(spacegroup, str): 

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

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

732 else: 

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

734 while True: 

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

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

737 _symbol = format_symbol(_symbol) 

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

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

740 _no = int(_no) 

741 

742 condition = ( 

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

744 and _setting == setting) 

745 or (isinstance(spacegroup, str) 

746 and compact_symbol == compact_spacegroup) and 

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

748 

749 if condition: 

750 namespace = SimpleNamespace() 

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

752 return namespace 

753 else: 

754 _skip_to_blank(f, spacegroup, setting) 

755 

756 

757def parse_sitesym_element(element): 

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

759 by the International Tables. 

760 

761 Examples 

762 -------- 

763 

764 >>> parse_sitesym_element("x") 

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

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

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

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

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

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

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

772 

773 

774 

775 Parameters 

776 ---------- 

777 

778 element: str 

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

780 

781 

782 Returns 

783 ------- 

784 

785 list[tuple[int, int]] 

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

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

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

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

790 

791 float 

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

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

794 

795 

796 """ 

797 element = element.lower() 

798 is_positive = True 

799 is_frac = False 

800 sng_trans = None 

801 fst_trans = [] 

802 snd_trans = [] 

803 rot = [] 

804 

805 for char in element: 

806 if char == "+": 

807 is_positive = True 

808 elif char == "-": 

809 is_positive = False 

810 elif char == "/": 

811 is_frac = True 

812 elif char in "xyz": 

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

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

815 if sng_trans is None: 

816 sng_trans = 1.0 if is_positive else -1.0 

817 if is_frac: 

818 snd_trans.append(char) 

819 else: 

820 fst_trans.append(char) 

821 

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

823 if is_frac: 

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

825 

826 return rot, trans 

827 

828 

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

830 force_positive_translation=False): 

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

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

833 

834 Parameters 

835 ---------- 

836 

837 sym: str 

838 Site symmetry in the form used by International Tables 

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

840 

841 out_rot: np.array 

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

843 

844 out_rot: np.array 

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

846 

847 sep: str 

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

849 

850 force_positive_translation: bool 

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

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

853 

854 

855 Returns 

856 ------- 

857 

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

859 

860 

861 """ 

862 out_rot[:] = 0.0 

863 out_trans[:] = 0.0 

864 

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

866 e_rot_list, e_trans = parse_sitesym_element(element) 

867 for rot_idx, rot_sgn in e_rot_list: 

868 out_rot[i][rot_idx] = rot_sgn 

869 out_trans[i] = \ 

870 (e_trans % 1.0) if force_positive_translation else e_trans 

871 

872 

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

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

875 International Tables and returns corresponding rotation and 

876 translation arrays. 

877 

878 Example: 

879 

880 >>> symlist = [ 

881 ... 'x,y,z', 

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

883 ... '-y,-x,-z', 

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

885 ... ] 

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

887 >>> rot 

888 array([[[ 1, 0, 0], 

889 [ 0, 1, 0], 

890 [ 0, 0, 1]], 

891 <BLANKLINE> 

892 [[ 0, -1, 0], 

893 [ 1, 0, 0], 

894 [ 0, 0, 1]], 

895 <BLANKLINE> 

896 [[ 0, -1, 0], 

897 [-1, 0, 0], 

898 [ 0, 0, -1]], 

899 <BLANKLINE> 

900 [[ 1, 0, 0], 

901 [ 0, 1, 0], 

902 [ 0, 0, -1]]]) 

903 >>> trans 

904 array([[ 0. , 0. , 0. ], 

905 [ 0.5 , 0.5 , 0. ], 

906 [ 0. , 0. , 0. ], 

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

908 """ 

909 

910 nsym = len(symlist) 

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

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

913 

914 for i, sym in enumerate(symlist): 

915 parse_sitesym_single( 

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

917 force_positive_translation=force_positive_translation) 

918 

919 return rot, trans 

920 

921 

922def spacegroup_from_data(no=None, 

923 symbol=None, 

924 setting=None, 

925 centrosymmetric=None, 

926 scaled_primitive_cell=None, 

927 reciprocal_cell=None, 

928 subtrans=None, 

929 sitesym=None, 

930 rotations=None, 

931 translations=None, 

932 datafile=None): 

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

934 useful when reading crystal data with its own spacegroup 

935 definitions.""" 

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

937 spg = Spacegroup(no, setting, datafile) 

938 elif symbol is not None: 

939 spg = Spacegroup(symbol, None, datafile) 

940 else: 

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

942 'or *symbol* must be given') 

943 if not isinstance(sitesym, list): 

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

945 

946 have_sym = False 

947 if centrosymmetric is not None: 

948 spg._centrosymmetric = bool(centrosymmetric) 

949 if scaled_primitive_cell is not None: 

950 spg._scaled_primitive_cell = np.array(scaled_primitive_cell) 

951 if reciprocal_cell is not None: 

952 spg._reciprocal_cell = np.array(reciprocal_cell) 

953 if subtrans is not None: 

954 spg._subtrans = np.atleast_2d(subtrans) 

955 if sitesym is not None: 

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

957 have_sym = True 

958 if rotations is not None: 

959 spg._rotations = np.atleast_3d(rotations) 

960 have_sym = True 

961 if translations is not None: 

962 spg._translations = np.atleast_2d(translations) 

963 have_sym = True 

964 if have_sym: 

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

966 raise SpacegroupValueError('inconsistent number of rotations and ' 

967 'translations') 

968 return spg 

969 

970 

971@deprecated( 

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

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

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

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

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

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

978) 

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

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

981 

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

983 

984 .. warning:: 

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

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

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

988 

989 .. deprecated:: 3.24.0 

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

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

992 

993 Parameters 

994 ---------- 

995 

996 atoms: Atoms object 

997 Types, positions and unit-cell. 

998 symprec: float 

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

1000 coordinates to find crystal symmetry. 

1001 

1002 The Spacegroup object is returned. 

1003 """ 

1004 

1005 # Example: 

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

1007 # when import fails) 

1008 # >>> from ase.build import bulk 

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

1010 # >>> sg = get_spacegroup(atoms) 

1011 # >>> sg 

1012 # Spacegroup(225, setting=1) 

1013 # >>> sg.no 

1014 # 225 

1015 

1016 import spglib 

1017 

1018 sg = spglib_new_errorhandling(spglib.get_spacegroup)( 

1019 (atoms.get_cell(), atoms.get_scaled_positions(), 

1020 atoms.get_atomic_numbers()), 

1021 symprec=symprec) 

1022 if sg is None: 

1023 raise RuntimeError('Spacegroup not found') 

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

1025 return Spacegroup(sg_no)