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

424 statements  

« prev     ^ index     » next       coverage.py v7.5.3, created at 2025-08-02 00:12 +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 

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 spacegroup : int | string | Spacegroup instance 

151 The space group number in International Tables of 

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

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

154 setting : 1 | 2 

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

156 determines Which of these should be used. 

157 datafile : None | string 

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

159 will be used. 

160 """ 

161 if isinstance(spacegroup, Spacegroup): 

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

163 setattr(self, k, v) 

164 return 

165 if not datafile: 

166 datafile = get_datafile() 

167 namespace = _read_datafile(spacegroup, setting, datafile) 

168 self._no = namespace._no 

169 self._symbol = namespace._symbol 

170 self._setting = namespace._setting 

171 self._centrosymmetric = namespace._centrosymmetric 

172 self._scaled_primitive_cell = namespace._scaled_primitive_cell 

173 self._reciprocal_cell = namespace._reciprocal_cell 

174 self._subtrans = namespace._subtrans 

175 self._rotations = namespace._rotations 

176 self._translations = namespace._translations 

177 

178 def __repr__(self): 

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

180 

181 def todict(self): 

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

183 

184 def __str__(self): 

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

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

187 retval = [] 

188 # no, symbol 

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

190 # setting 

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

192 # centrosymmetric 

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

194 # primitive vectors 

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

196 for i in range(3): 

197 retval.append(' ') 

198 for j in range(3): 

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

200 retval.append('\n') 

201 # primitive reciprocal vectors 

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

203 for i in range(3): 

204 retval.append(' ') 

205 for j in range(3): 

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

207 retval.append('\n') 

208 # sublattice 

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

210 for i in range(self.nsubtrans): 

211 retval.append(' ') 

212 for j in range(3): 

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

214 retval.append('\n') 

215 # symmetry operations 

216 nrot = len(self.rotations) 

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

218 for i in range(nrot): 

219 retval.append(' ') 

220 for j in range(3): 

221 retval.append(' ') 

222 for k in range(3): 

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

224 retval.append(' ') 

225 for j in range(3): 

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

227 retval.append('\n') 

228 retval.append('\n') 

229 return ''.join(retval) 

230 

231 def __eq__(self, other): 

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

233 

234 def __ne__(self, other): 

235 return not self.__eq__(other) 

236 

237 def __lt__(self, other): 

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

239 and self.setting < other.setting) 

240 

241 def __index__(self): 

242 return self.no 

243 

244 __int__ = __index__ 

245 

246 def get_symop(self): 

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

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

249 tuples.""" 

250 symop = [] 

251 parities = [1] 

252 if self.centrosymmetric: 

253 parities.append(-1) 

254 for parity in parities: 

255 for subtrans in self.subtrans: 

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

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

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

259 return symop 

260 

261 def get_op(self): 

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

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

264 two ndarrays.""" 

265 if self.centrosymmetric: 

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

267 (self.nsubtrans, 1, 1)) 

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

269 (self.nsubtrans, 1)) 

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

271 trans = np.mod(trans, 1) 

272 else: 

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

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

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

276 trans = np.mod(trans, 1) 

277 return rot, trans 

278 

279 def get_rotations(self): 

280 """Return all rotations, including inversions for 

281 centrosymmetric crystals.""" 

282 if self.centrosymmetric: 

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

284 else: 

285 return self.rotations 

286 

287 def equivalent_reflections(self, hkl): 

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

289 in hkl. 

290 

291 Example: 

292 

293 >>> from ase.spacegroup import Spacegroup 

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

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

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

297 [ 0, -2, 0], 

298 [-2, 0, 0], 

299 [ 2, 0, 0], 

300 [ 0, 2, 0], 

301 [ 0, 0, 2]]) 

302 """ 

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

304 rot = self.get_rotations() 

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

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

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

308 ind = np.lexsort(refl.T) 

309 refl = refl[ind] 

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

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

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

313 

314 def equivalent_lattice_points(self, uvw): 

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

316 in `uvw` with respect to rotations only. 

317 

318 Only equivalent lattice points that conserves the distance to 

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

320 space version of the equivalent_reflections() method). 

321 

322 Example: 

323 

324 >>> from ase.spacegroup import Spacegroup 

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

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

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

328 [ 0, -2, 0], 

329 [-2, 0, 0], 

330 [ 2, 0, 0], 

331 [ 0, 2, 0], 

332 [ 0, 0, 2]]) 

333 

334 """ 

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

336 rot = self.get_rotations() 

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

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

339 ind = np.lexsort(directions.T) 

340 directions = directions[ind] 

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

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

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

344 

345 def symmetry_normalised_reflections(self, hkl): 

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

347 corresponding symmetry-equivalent reflections of lowest 

348 indices. 

349 

350 Example: 

351 

352 >>> from ase.spacegroup import Spacegroup 

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

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

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

356 [ 0, 0, -2]]) 

357 """ 

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

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

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

361 for i, g in enumerate(hkl): 

362 gsym = np.dot(R, g) 

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

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

365 return normalised 

366 

367 def unique_reflections(self, hkl): 

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

369 reflections. 

370 

371 Example: 

372 

373 >>> from ase.spacegroup import Spacegroup 

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

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

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

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

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

379 array([[2, 0, 0], 

380 [2, 2, 0]]) 

381 """ 

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

383 hklnorm = self.symmetry_normalised_reflections(hkl) 

384 perm = np.lexsort(hklnorm.T) 

385 iperm = perm.argsort() 

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

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

388 imask = mask[iperm] 

389 return hkl[imask] 

390 

391 def equivalent_sites(self, 

392 scaled_positions, 

393 onduplicates='error', 

394 symprec=1e-3, 

395 occupancies=None): 

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

397 

398 Parameters: 

399 

400 scaled_positions: list | array 

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

402 

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

404 List of occupancies corresponding to the respective sites. 

405 

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

407 Action if `scaled_positions` contain symmetry-equivalent 

408 positions of full occupancy: 

409 

410 'keep' 

411 ignore additional symmetry-equivalent positions 

412 'replace' 

413 replace 

414 'warn' 

415 like 'keep', but issue an UserWarning 

416 'error' 

417 raises a SpacegroupValueError 

418 

419 symprec: float 

420 Minimum "distance" betweed two sites in scaled coordinates 

421 before they are counted as the same site. 

422 

423 Returns: 

424 

425 sites: array 

426 A NumPy array of equivalent sites. 

427 kinds: list 

428 A list of integer indices specifying which input site is 

429 equivalent to the corresponding returned site. 

430 

431 Example: 

432 

433 >>> from ase.spacegroup import Spacegroup 

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

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

436 >>> sites 

437 array([[ 0. , 0. , 0. ], 

438 [ 0. , 0.5, 0.5], 

439 [ 0.5, 0. , 0.5], 

440 [ 0.5, 0.5, 0. ], 

441 [ 0.5, 0. , 0. ], 

442 [ 0. , 0.5, 0. ], 

443 [ 0. , 0. , 0.5], 

444 [ 0.5, 0.5, 0.5]]) 

445 >>> kinds 

446 [0, 0, 0, 0, 1, 1, 1, 1] 

447 """ 

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

449 raise SpacegroupValueError( 

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

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

452 ) 

453 

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

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

456 rotations = np.array(rotations) 

457 translations = np.array(translations) 

458 

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

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

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

462 orbit = [candidates[0]] 

463 for member in candidates[1:]: 

464 diff = member - orbit 

465 diff -= np.rint(diff) 

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

467 orbit.append(member) 

468 return np.array(orbit) 

469 

470 orbits = [] 

471 for kind, pos in enumerate(scaled): 

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

473 diff = pos - positions0 

474 diff -= np.rint(diff) 

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

476 if onduplicates == 'keep': 

477 pass 

478 elif onduplicates == 'replace': 

479 orbits[i] = (kind, positions0) 

480 elif onduplicates == 'warn': 

481 warnings.warn( 

482 'scaled_positions %d and %d are equivalent' % 

483 (kind0, kind)) 

484 elif onduplicates == 'error': 

485 raise SpacegroupValueError( 

486 'scaled_positions %d and %d are equivalent' % 

487 (kind0, kind)) 

488 break 

489 else: 

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

491 

492 kinds = [] 

493 sites = [] 

494 for kind, orbit in orbits: 

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

496 sites.append(orbit) 

497 

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

499 

500 def symmetry_normalised_sites(self, 

501 scaled_positions, 

502 map_to_unitcell=True): 

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

504 containing the corresponding symmetry-equivalent sites of 

505 lowest indices. 

506 

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

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

509 included as symmetry operator. 

510 

511 Example: 

512 

513 >>> from ase.spacegroup import Spacegroup 

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

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

516 array([[ 0., 0., 0.], 

517 [ 0., 0., 0.]]) 

518 """ 

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

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

521 rot, trans = self.get_op() 

522 for i, pos in enumerate(scaled): 

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

524 if map_to_unitcell: 

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

526 sympos %= 1.0 

527 sympos %= 1.0 

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

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

530 return normalised 

531 

532 def unique_sites(self, 

533 scaled_positions, 

534 symprec=1e-3, 

535 output_mask=False, 

536 map_to_unitcell=True): 

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

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

539 array masking the subset is also returned. 

540 

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

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

543 

544 Example: 

545 

546 >>> from ase.spacegroup import Spacegroup 

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

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

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

550 ... [1.0, 0.0, 0.0], 

551 ... [0.5, 0.0, 0.0]]) 

552 array([[ 0. , 0. , 0. ], 

553 [ 0.5, 0. , 0. ]]) 

554 """ 

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

556 symnorm = self.symmetry_normalised_sites(scaled, map_to_unitcell) 

557 perm = np.lexsort(symnorm.T) 

558 iperm = perm.argsort() 

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

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

561 imask = mask[iperm] 

562 if output_mask: 

563 return scaled[imask], imask 

564 else: 

565 return scaled[imask] 

566 

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

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

569 tagging all equivalent atoms with the same index. 

570 

571 Example: 

572 

573 >>> from ase.spacegroup import Spacegroup 

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

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

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

577 ... [1.0, 0.0, 0.0], 

578 ... [0.5, 0.0, 0.0]]) 

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

580 """ 

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

582 scaled %= 1.0 

583 scaled %= 1.0 

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

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

586 rot, trans = self.get_op() 

587 i = 0 

588 while mask.any(): 

589 pos = scaled[mask][0] 

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

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

592 sympos %= 1.0 

593 sympos %= 1.0 

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

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

596 axis=2), 

597 axis=0) 

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

599 tags[m] = i 

600 mask &= ~m 

601 i += 1 

602 return tags 

603 

604 

605def get_datafile(): 

606 """Return default path to datafile.""" 

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

608 

609 

610def format_symbol(symbol): 

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

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

613 removing dublicated spaces.""" 

614 fixed = [] 

615 s = symbol.strip() 

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

617 for c in s: 

618 if c.isalpha(): 

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

620 fixed.append(c) 

621 else: 

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

623 elif c.isspace(): 

624 fixed.append(' ') 

625 elif c.isdigit(): 

626 fixed.append(c) 

627 elif c == '-': 

628 fixed.append(' ' + c) 

629 elif c == '/': 

630 fixed.append(c) 

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

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

633 

634 

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

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

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

638# instance is created. 

639 

640 

641def _skip_to_blank(f, spacegroup, setting): 

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

643 while True: 

644 line = f.readline() 

645 if not line: 

646 raise SpacegroupNotFoundError( 

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

648 'found in data base') 

649 if not line.strip(): 

650 break 

651 

652 

653def _skip_to_nonblank(f, spacegroup, setting): 

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

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

656 while True: 

657 line1 = f.readline() 

658 if not line1: 

659 raise SpacegroupNotFoundError( 

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

661 (spacegroup, setting)) 

662 line1.strip() 

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

664 line2 = f.readline() 

665 break 

666 return line1, line2 

667 

668 

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

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

671 

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

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

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

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

676 

677 spg._no = no 

678 spg._symbol = symbol.strip() 

679 spg._setting = setting 

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

681 # primitive vectors 

682 f.readline() 

683 spg._scaled_primitive_cell = np.array( 

684 [ 

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

686 for _ in range(3) 

687 ], 

688 dtype=float, 

689 ) 

690 # primitive reciprocal vectors 

691 f.readline() 

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

693 for i in range(3)], 

694 dtype=int) 

695 # subtranslations 

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

697 spg._subtrans = np.array( 

698 [ 

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

700 for _ in range(nsubtrans) 

701 ], 

702 dtype=float, 

703 ) 

704 # symmetry operations 

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

706 symop = np.array( 

707 [ 

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

709 for _ in range(nsym) 

710 ], 

711 dtype=float, 

712 ) 

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

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

715 

716 

717@lru_cache 

718def _read_datafile(spacegroup, setting, datafile): 

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

720 return _read_f(spacegroup, setting, fd) 

721 

722 

723def _read_f(spacegroup, setting, f): 

724 if isinstance(spacegroup, int): 

725 pass 

726 elif isinstance(spacegroup, str): 

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

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

729 else: 

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

731 while True: 

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

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

734 _symbol = format_symbol(_symbol) 

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

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

737 _no = int(_no) 

738 

739 condition = ( 

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

741 and _setting == setting) 

742 or (isinstance(spacegroup, str) 

743 and compact_symbol == compact_spacegroup) and 

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

745 

746 if condition: 

747 namespace = SimpleNamespace() 

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

749 return namespace 

750 else: 

751 _skip_to_blank(f, spacegroup, setting) 

752 

753 

754def parse_sitesym_element(element): 

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

756 by the International Tables. 

757 

758 Examples: 

759 

760 >>> parse_sitesym_element("x") 

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

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

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

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

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

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

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

768 

769 

770 

771 Parameters 

772 ---------- 

773 

774 element: str 

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

776 

777 

778 Returns 

779 ------- 

780 

781 list[tuple[int, int]] 

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

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

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

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

786 

787 float 

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

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

790 

791 

792 """ 

793 element = element.lower() 

794 is_positive = True 

795 is_frac = False 

796 sng_trans = None 

797 fst_trans = [] 

798 snd_trans = [] 

799 rot = [] 

800 

801 for char in element: 

802 if char == "+": 

803 is_positive = True 

804 elif char == "-": 

805 is_positive = False 

806 elif char == "/": 

807 is_frac = True 

808 elif char in "xyz": 

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

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

811 if sng_trans is None: 

812 sng_trans = 1.0 if is_positive else -1.0 

813 if is_frac: 

814 snd_trans.append(char) 

815 else: 

816 fst_trans.append(char) 

817 

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

819 if is_frac: 

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

821 

822 return rot, trans 

823 

824 

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

826 force_positive_translation=False): 

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

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

829 

830 Parameters 

831 ---------- 

832 

833 sym: str 

834 Site symmetry in the form used by International Tables 

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

836 

837 out_rot: np.array 

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

839 

840 out_rot: np.array 

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

842 

843 sep: str 

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

845 

846 force_positive_translation: bool 

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

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

849 

850 

851 Returns 

852 ------- 

853 

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

855 

856 

857 """ 

858 out_rot[:] = 0.0 

859 out_trans[:] = 0.0 

860 

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

862 e_rot_list, e_trans = parse_sitesym_element(element) 

863 for rot_idx, rot_sgn in e_rot_list: 

864 out_rot[i][rot_idx] = rot_sgn 

865 out_trans[i] = \ 

866 (e_trans % 1.0) if force_positive_translation else e_trans 

867 

868 

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

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

871 International Tables and returns corresponding rotation and 

872 translation arrays. 

873 

874 Example: 

875 

876 >>> symlist = [ 

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

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

879 ... '-y,-x,-z', 

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

881 ... ] 

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

883 >>> rot 

884 array([[[ 1, 0, 0], 

885 [ 0, 1, 0], 

886 [ 0, 0, 1]], 

887 <BLANKLINE> 

888 [[ 0, -1, 0], 

889 [ 1, 0, 0], 

890 [ 0, 0, 1]], 

891 <BLANKLINE> 

892 [[ 0, -1, 0], 

893 [-1, 0, 0], 

894 [ 0, 0, -1]], 

895 <BLANKLINE> 

896 [[ 1, 0, 0], 

897 [ 0, 1, 0], 

898 [ 0, 0, -1]]]) 

899 >>> trans 

900 array([[ 0. , 0. , 0. ], 

901 [ 0.5 , 0.5 , 0. ], 

902 [ 0. , 0. , 0. ], 

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

904 """ 

905 

906 nsym = len(symlist) 

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

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

909 

910 for i, sym in enumerate(symlist): 

911 parse_sitesym_single( 

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

913 force_positive_translation=force_positive_translation) 

914 

915 return rot, trans 

916 

917 

918def spacegroup_from_data(no=None, 

919 symbol=None, 

920 setting=None, 

921 centrosymmetric=None, 

922 scaled_primitive_cell=None, 

923 reciprocal_cell=None, 

924 subtrans=None, 

925 sitesym=None, 

926 rotations=None, 

927 translations=None, 

928 datafile=None): 

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

930 useful when reading crystal data with its own spacegroup 

931 definitions.""" 

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

933 spg = Spacegroup(no, setting, datafile) 

934 elif symbol is not None: 

935 spg = Spacegroup(symbol, None, datafile) 

936 else: 

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

938 'or *symbol* must be given') 

939 if not isinstance(sitesym, list): 

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

941 

942 have_sym = False 

943 if centrosymmetric is not None: 

944 spg._centrosymmetric = bool(centrosymmetric) 

945 if scaled_primitive_cell is not None: 

946 spg._scaled_primitive_cell = np.array(scaled_primitive_cell) 

947 if reciprocal_cell is not None: 

948 spg._reciprocal_cell = np.array(reciprocal_cell) 

949 if subtrans is not None: 

950 spg._subtrans = np.atleast_2d(subtrans) 

951 if sitesym is not None: 

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

953 have_sym = True 

954 if rotations is not None: 

955 spg._rotations = np.atleast_3d(rotations) 

956 have_sym = True 

957 if translations is not None: 

958 spg._translations = np.atleast_2d(translations) 

959 have_sym = True 

960 if have_sym: 

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

962 raise SpacegroupValueError('inconsistent number of rotations and ' 

963 'translations') 

964 return spg 

965 

966 

967@deprecated( 

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

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

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

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

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

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

974) 

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

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

977 

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

979 

980 .. warning:: 

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

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

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

984 

985 .. deprecated:: 3.24.0 

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

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

988 

989 Parameters: 

990 

991 atoms: Atoms object 

992 Types, positions and unit-cell. 

993 symprec: float 

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

995 coordinates to find crystal symmetry. 

996 

997 The Spacegroup object is returned. 

998 """ 

999 

1000 # Example: 

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

1002 # when import fails) 

1003 # >>> from ase.build import bulk 

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

1005 # >>> sg = get_spacegroup(atoms) 

1006 # >>> sg 

1007 # Spacegroup(225, setting=1) 

1008 # >>> sg.no 

1009 # 225 

1010 

1011 import spglib 

1012 

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

1014 atoms.get_atomic_numbers()), 

1015 symprec=symprec) 

1016 if sg is None: 

1017 raise RuntimeError('Spacegroup not found') 

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

1019 return Spacegroup(sg_no)