Coverage for /builds/debichem-team/python-ase/ase/spacegroup/spacegroup.py : 92.15%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Copyright (C) 2010, Jesper Friis
2# (see accompanying license files for details).
3"""Definition of the Spacegroup class.
5This module only depends on NumPy and the space group database.
6"""
8import os
9import warnings
10from functools import total_ordering
11from typing import Union
13import numpy as np
15__all__ = ['Spacegroup']
18class SpacegroupError(Exception):
19 """Base exception for the spacegroup module."""
20 pass
23class SpacegroupNotFoundError(SpacegroupError):
24 """Raised when given space group cannot be found in data base."""
25 pass
28class SpacegroupValueError(SpacegroupError):
29 """Raised when arguments have invalid value."""
30 pass
33# Type alias
34_SPACEGROUP = Union[int, str, 'Spacegroup']
37@total_ordering
38class Spacegroup:
39 """A space group class.
41 The instances of Spacegroup describes the symmetry operations for
42 the given space group.
44 Example:
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 no = property(
63 lambda self: self._no,
64 doc='Space group number in International Tables of Crystallography.')
65 symbol = property(
66 lambda self: self._symbol,
67 doc='Hermann-Mauguin (or international) symbol for the space group.')
68 setting = property(lambda self: self._setting,
69 doc='Space group setting. Either one or two.')
70 lattice = property(lambda self: self._symbol[0],
71 doc="""Lattice type:
73 P primitive
74 I body centering, h+k+l=2n
75 F face centering, h,k,l all odd or even
76 A,B,C single face centering, k+l=2n, h+l=2n, h+k=2n
77 R rhombohedral centering, -h+k+l=3n (obverse); h-k+l=3n (reverse)
78 """)
79 centrosymmetric = property(lambda self: self._centrosymmetric,
80 doc='Whether a center of symmetry exists.')
81 scaled_primitive_cell = property(
82 lambda self: self._scaled_primitive_cell,
83 doc='Primitive cell in scaled coordinates as a matrix with the '
84 'primitive vectors along the rows.')
85 reciprocal_cell = property(
86 lambda self: self._reciprocal_cell,
87 doc='Tree Miller indices that span all kinematically non-forbidden '
88 'reflections as a matrix with the Miller indices along the rows.')
89 nsubtrans = property(lambda self: len(self._subtrans),
90 doc='Number of cell-subtranslation vectors.')
92 def _get_nsymop(self):
93 """Returns total number of symmetry operations."""
94 if self.centrosymmetric:
95 return 2 * len(self._rotations) * len(self._subtrans)
96 else:
97 return len(self._rotations) * len(self._subtrans)
99 nsymop = property(_get_nsymop, doc='Total number of symmetry operations.')
100 subtrans = property(
101 lambda self: self._subtrans,
102 doc='Translations vectors belonging to cell-sub-translations.')
103 rotations = property(
104 lambda self: self._rotations,
105 doc='Symmetry rotation matrices. The invertions are not included '
106 'for centrosymmetrical crystals.')
107 translations = property(
108 lambda self: self._translations,
109 doc='Symmetry translations. The invertions are not included '
110 'for centrosymmetrical crystals.')
112 def __init__(self, spacegroup: _SPACEGROUP, setting=1, datafile=None):
113 """Returns a new Spacegroup instance.
115 Parameters:
117 spacegroup : int | string | Spacegroup instance
118 The space group number in International Tables of
119 Crystallography or its Hermann-Mauguin symbol. E.g.
120 spacegroup=225 and spacegroup='F m -3 m' are equivalent.
121 setting : 1 | 2
122 Some space groups have more than one setting. `setting`
123 determines Which of these should be used.
124 datafile : None | string
125 Path to database file. If `None`, the the default database
126 will be used.
127 """
128 if isinstance(spacegroup, Spacegroup):
129 for k, v in spacegroup.__dict__.items():
130 setattr(self, k, v)
131 return
132 if not datafile:
133 datafile = get_datafile()
134 with open(datafile, 'r') as fd:
135 _read_datafile(self, spacegroup, setting, fd)
137 def __repr__(self):
138 return 'Spacegroup(%d, setting=%d)' % (self.no, self.setting)
140 def todict(self):
141 return {'number': self.no, 'setting': self.setting}
143 def __str__(self):
144 """Return a string representation of the space group data in
145 the same format as found the database."""
146 retval = []
147 # no, symbol
148 retval.append('%-3d %s\n' % (self.no, self.symbol))
149 # setting
150 retval.append(' setting %d\n' % (self.setting))
151 # centrosymmetric
152 retval.append(' centrosymmetric %d\n' % (self.centrosymmetric))
153 # primitive vectors
154 retval.append(' primitive vectors\n')
155 for i in range(3):
156 retval.append(' ')
157 for j in range(3):
158 retval.append(' %13.10f' % (self.scaled_primitive_cell[i, j]))
159 retval.append('\n')
160 # primitive reciprocal vectors
161 retval.append(' reciprocal vectors\n')
162 for i in range(3):
163 retval.append(' ')
164 for j in range(3):
165 retval.append(' %3d' % (self.reciprocal_cell[i, j]))
166 retval.append('\n')
167 # sublattice
168 retval.append(' %d subtranslations\n' % self.nsubtrans)
169 for i in range(self.nsubtrans):
170 retval.append(' ')
171 for j in range(3):
172 retval.append(' %13.10f' % (self.subtrans[i, j]))
173 retval.append('\n')
174 # symmetry operations
175 nrot = len(self.rotations)
176 retval.append(' %d symmetry operations (rot+trans)\n' % nrot)
177 for i in range(nrot):
178 retval.append(' ')
179 for j in range(3):
180 retval.append(' ')
181 for k in range(3):
182 retval.append(' %2d' % (self.rotations[i, j, k]))
183 retval.append(' ')
184 for j in range(3):
185 retval.append(' %13.10f' % self.translations[i, j])
186 retval.append('\n')
187 retval.append('\n')
188 return ''.join(retval)
190 def __eq__(self, other):
191 return self.no == other.no and self.setting == other.setting
193 def __ne__(self, other):
194 return not self.__eq__(other)
196 def __lt__(self, other):
197 return self.no < other.no or (self.no == other.no
198 and self.setting < other.setting)
200 def __index__(self):
201 return self.no
203 __int__ = __index__
205 def get_symop(self):
206 """Returns all symmetry operations (including inversions and
207 subtranslations) as a sequence of (rotation, translation)
208 tuples."""
209 symop = []
210 parities = [1]
211 if self.centrosymmetric:
212 parities.append(-1)
213 for parity in parities:
214 for subtrans in self.subtrans:
215 for rot, trans in zip(self.rotations, self.translations):
216 newtrans = np.mod(trans + subtrans, 1)
217 symop.append((parity * rot, newtrans))
218 return symop
220 def get_op(self):
221 """Returns all symmetry operations (including inversions and
222 subtranslations), but unlike get_symop(), they are returned as
223 two ndarrays."""
224 if self.centrosymmetric:
225 rot = np.tile(np.vstack((self.rotations, -self.rotations)),
226 (self.nsubtrans, 1, 1))
227 trans = np.tile(np.vstack((self.translations, -self.translations)),
228 (self.nsubtrans, 1))
229 trans += np.repeat(self.subtrans, 2 * len(self.rotations), axis=0)
230 trans = np.mod(trans, 1)
231 else:
232 rot = np.tile(self.rotations, (self.nsubtrans, 1, 1))
233 trans = np.tile(self.translations, (self.nsubtrans, 1))
234 trans += np.repeat(self.subtrans, len(self.rotations), axis=0)
235 trans = np.mod(trans, 1)
236 return rot, trans
238 def get_rotations(self):
239 """Return all rotations, including inversions for
240 centrosymmetric crystals."""
241 if self.centrosymmetric:
242 return np.vstack((self.rotations, -self.rotations))
243 else:
244 return self.rotations
246 def equivalent_reflections(self, hkl):
247 """Return all equivalent reflections to the list of Miller indices
248 in hkl.
250 Example:
252 >>> from ase.spacegroup import Spacegroup
253 >>> sg = Spacegroup(225) # fcc
254 >>> sg.equivalent_reflections([[0, 0, 2]])
255 array([[ 0, 0, -2],
256 [ 0, -2, 0],
257 [-2, 0, 0],
258 [ 2, 0, 0],
259 [ 0, 2, 0],
260 [ 0, 0, 2]])
261 """
262 hkl = np.array(hkl, dtype='int', ndmin=2)
263 rot = self.get_rotations()
264 n, nrot = len(hkl), len(rot)
265 R = rot.transpose(0, 2, 1).reshape((3 * nrot, 3)).T
266 refl = np.dot(hkl, R).reshape((n * nrot, 3))
267 ind = np.lexsort(refl.T)
268 refl = refl[ind]
269 diff = np.diff(refl, axis=0)
270 mask = np.any(diff, axis=1)
271 return np.vstack((refl[:-1][mask], refl[-1, :]))
273 def equivalent_lattice_points(self, uvw):
274 """Return all lattice points equivalent to any of the lattice points
275 in `uvw` with respect to rotations only.
277 Only equivalent lattice points that conserves the distance to
278 origo are included in the output (making this a kind of real
279 space version of the equivalent_reflections() method).
281 Example:
283 >>> from ase.spacegroup import Spacegroup
284 >>> sg = Spacegroup(225) # fcc
285 >>> sg.equivalent_lattice_points([[0, 0, 2]])
286 array([[ 0, 0, -2],
287 [ 0, -2, 0],
288 [-2, 0, 0],
289 [ 2, 0, 0],
290 [ 0, 2, 0],
291 [ 0, 0, 2]])
293 """
294 uvw = np.array(uvw, ndmin=2)
295 rot = self.get_rotations()
296 n, nrot = len(uvw), len(rot)
297 directions = np.dot(uvw, rot).reshape((n * nrot, 3))
298 ind = np.lexsort(directions.T)
299 directions = directions[ind]
300 diff = np.diff(directions, axis=0)
301 mask = np.any(diff, axis=1)
302 return np.vstack((directions[:-1][mask], directions[-1:]))
304 def symmetry_normalised_reflections(self, hkl):
305 """Returns an array of same size as *hkl*, containing the
306 corresponding symmetry-equivalent reflections of lowest
307 indices.
309 Example:
311 >>> from ase.spacegroup import Spacegroup
312 >>> sg = Spacegroup(225) # fcc
313 >>> sg.symmetry_normalised_reflections([[2, 0, 0], [0, 2, 0]])
314 array([[ 0, 0, -2],
315 [ 0, 0, -2]])
316 """
317 hkl = np.array(hkl, dtype=int, ndmin=2)
318 normalised = np.empty(hkl.shape, int)
319 R = self.get_rotations().transpose(0, 2, 1)
320 for i, g in enumerate(hkl):
321 gsym = np.dot(R, g)
322 j = np.lexsort(gsym.T)[0]
323 normalised[i, :] = gsym[j]
324 return normalised
326 def unique_reflections(self, hkl):
327 """Returns a subset *hkl* containing only the symmetry-unique
328 reflections.
330 Example:
332 >>> from ase.spacegroup import Spacegroup
333 >>> sg = Spacegroup(225) # fcc
334 >>> sg.unique_reflections([[ 2, 0, 0],
335 ... [ 0, -2, 0],
336 ... [ 2, 2, 0],
337 ... [ 0, -2, -2]])
338 array([[2, 0, 0],
339 [2, 2, 0]])
340 """
341 hkl = np.array(hkl, dtype=int, ndmin=2)
342 hklnorm = self.symmetry_normalised_reflections(hkl)
343 perm = np.lexsort(hklnorm.T)
344 iperm = perm.argsort()
345 xmask = np.abs(np.diff(hklnorm[perm], axis=0)).any(axis=1)
346 mask = np.concatenate(([True], xmask))
347 imask = mask[iperm]
348 return hkl[imask]
350 def equivalent_sites(self,
351 scaled_positions,
352 onduplicates='error',
353 symprec=1e-3,
354 occupancies=None):
355 """Returns the scaled positions and all their equivalent sites.
357 Parameters:
359 scaled_positions: list | array
360 List of non-equivalent sites given in unit cell coordinates.
362 occupancies: list | array, optional (default=None)
363 List of occupancies corresponding to the respective sites.
365 onduplicates : 'keep' | 'replace' | 'warn' | 'error'
366 Action if `scaled_positions` contain symmetry-equivalent
367 positions of full occupancy:
369 'keep'
370 ignore additional symmetry-equivalent positions
371 'replace'
372 replace
373 'warn'
374 like 'keep', but issue an UserWarning
375 'error'
376 raises a SpacegroupValueError
378 symprec: float
379 Minimum "distance" betweed two sites in scaled coordinates
380 before they are counted as the same site.
382 Returns:
384 sites: array
385 A NumPy array of equivalent sites.
386 kinds: list
387 A list of integer indices specifying which input site is
388 equivalent to the corresponding returned site.
390 Example:
392 >>> from ase.spacegroup import Spacegroup
393 >>> sg = Spacegroup(225) # fcc
394 >>> sites, kinds = sg.equivalent_sites([[0, 0, 0], [0.5, 0.0, 0.0]])
395 >>> sites
396 array([[ 0. , 0. , 0. ],
397 [ 0. , 0.5, 0.5],
398 [ 0.5, 0. , 0.5],
399 [ 0.5, 0.5, 0. ],
400 [ 0.5, 0. , 0. ],
401 [ 0. , 0.5, 0. ],
402 [ 0. , 0. , 0.5],
403 [ 0.5, 0.5, 0.5]])
404 >>> kinds
405 [0, 0, 0, 0, 1, 1, 1, 1]
406 """
407 kinds = []
408 sites = []
410 scaled = np.array(scaled_positions, ndmin=2)
412 for kind, pos in enumerate(scaled):
413 for rot, trans in self.get_symop():
414 site = np.mod(np.dot(rot, pos) + trans, 1.)
415 if not sites:
416 sites.append(site)
417 kinds.append(kind)
418 continue
419 t = site - sites
420 mask = np.all(
421 (abs(t) < symprec) | (abs(abs(t) - 1.0) < symprec), axis=1)
422 if np.any(mask):
423 inds = np.argwhere(mask).flatten()
424 for ind in inds:
425 # then we would just add the same thing again -> skip
426 if kinds[ind] == kind:
427 pass
428 elif onduplicates == 'keep':
429 pass
430 elif onduplicates == 'replace':
431 kinds[ind] = kind
432 elif onduplicates == 'warn':
433 warnings.warn('scaled_positions %d and %d '
434 'are equivalent' %
435 (kinds[ind], kind))
436 elif onduplicates == 'error':
437 raise SpacegroupValueError(
438 'scaled_positions %d and %d are equivalent' %
439 (kinds[ind], kind))
440 else:
441 raise SpacegroupValueError(
442 'Argument "onduplicates" must be one of: '
443 '"keep", "replace", "warn" or "error".')
444 else:
445 sites.append(site)
446 kinds.append(kind)
448 return np.array(sites), kinds
450 def symmetry_normalised_sites(self,
451 scaled_positions,
452 map_to_unitcell=True):
453 """Returns an array of same size as *scaled_positions*,
454 containing the corresponding symmetry-equivalent sites of
455 lowest indices.
457 If *map_to_unitcell* is true, the returned positions are all
458 mapped into the unit cell, i.e. lattice translations are
459 included as symmetry operator.
461 Example:
463 >>> from ase.spacegroup import Spacegroup
464 >>> sg = Spacegroup(225) # fcc
465 >>> sg.symmetry_normalised_sites([[0.0, 0.5, 0.5], [1.0, 1.0, 0.0]])
466 array([[ 0., 0., 0.],
467 [ 0., 0., 0.]])
468 """
469 scaled = np.array(scaled_positions, ndmin=2)
470 normalised = np.empty(scaled.shape, float)
471 rot, trans = self.get_op()
472 for i, pos in enumerate(scaled):
473 sympos = np.dot(rot, pos) + trans
474 if map_to_unitcell:
475 # Must be done twice, see the scaled_positions.py test
476 sympos %= 1.0
477 sympos %= 1.0
478 j = np.lexsort(sympos.T)[0]
479 normalised[i, :] = sympos[j]
480 return normalised
482 def unique_sites(self,
483 scaled_positions,
484 symprec=1e-3,
485 output_mask=False,
486 map_to_unitcell=True):
487 """Returns a subset of *scaled_positions* containing only the
488 symmetry-unique positions. If *output_mask* is True, a boolean
489 array masking the subset is also returned.
491 If *map_to_unitcell* is true, all sites are first mapped into
492 the unit cell making e.g. [0, 0, 0] and [1, 0, 0] equivalent.
494 Example:
496 >>> from ase.spacegroup import Spacegroup
497 >>> sg = Spacegroup(225) # fcc
498 >>> sg.unique_sites([[0.0, 0.0, 0.0],
499 ... [0.5, 0.5, 0.0],
500 ... [1.0, 0.0, 0.0],
501 ... [0.5, 0.0, 0.0]])
502 array([[ 0. , 0. , 0. ],
503 [ 0.5, 0. , 0. ]])
504 """
505 scaled = np.array(scaled_positions, ndmin=2)
506 symnorm = self.symmetry_normalised_sites(scaled, map_to_unitcell)
507 perm = np.lexsort(symnorm.T)
508 iperm = perm.argsort()
509 xmask = np.abs(np.diff(symnorm[perm], axis=0)).max(axis=1) > symprec
510 mask = np.concatenate(([True], xmask))
511 imask = mask[iperm]
512 if output_mask:
513 return scaled[imask], imask
514 else:
515 return scaled[imask]
517 def tag_sites(self, scaled_positions, symprec=1e-3):
518 """Returns an integer array of the same length as *scaled_positions*,
519 tagging all equivalent atoms with the same index.
521 Example:
523 >>> from ase.spacegroup import Spacegroup
524 >>> sg = Spacegroup(225) # fcc
525 >>> sg.tag_sites([[0.0, 0.0, 0.0],
526 ... [0.5, 0.5, 0.0],
527 ... [1.0, 0.0, 0.0],
528 ... [0.5, 0.0, 0.0]])
529 array([0, 0, 0, 1])
530 """
531 scaled = np.array(scaled_positions, ndmin=2)
532 scaled %= 1.0
533 scaled %= 1.0
534 tags = -np.ones((len(scaled), ), dtype=int)
535 mask = np.ones((len(scaled), ), dtype=bool)
536 rot, trans = self.get_op()
537 i = 0
538 while mask.any():
539 pos = scaled[mask][0]
540 sympos = np.dot(rot, pos) + trans
541 # Must be done twice, see the scaled_positions.py test
542 sympos %= 1.0
543 sympos %= 1.0
544 m = ~np.all(np.any(np.abs(scaled[np.newaxis, :, :] -
545 sympos[:, np.newaxis, :]) > symprec,
546 axis=2),
547 axis=0)
548 assert not np.any((~mask) & m)
549 tags[m] = i
550 mask &= ~m
551 i += 1
552 return tags
555def get_datafile():
556 """Return default path to datafile."""
557 return os.path.join(os.path.dirname(__file__), 'spacegroup.dat')
560def format_symbol(symbol):
561 """Returns well formatted Hermann-Mauguin symbol as extected by
562 the database, by correcting the case and adding missing or
563 removing dublicated spaces."""
564 fixed = []
565 s = symbol.strip()
566 s = s[0].upper() + s[1:].lower()
567 for c in s:
568 if c.isalpha():
569 if len(fixed) and fixed[-1] == '/':
570 fixed.append(c)
571 else:
572 fixed.append(' ' + c + ' ')
573 elif c.isspace():
574 fixed.append(' ')
575 elif c.isdigit():
576 fixed.append(c)
577 elif c == '-':
578 fixed.append(' ' + c)
579 elif c == '/':
580 fixed.append(c)
581 s = ''.join(fixed).strip()
582 return ' '.join(s.split())
585# Functions for parsing the database. They are moved outside the
586# Spacegroup class in order to make it easier to later implement
587# caching to avoid reading the database each time a new Spacegroup
588# instance is created.
591def _skip_to_blank(f, spacegroup, setting):
592 """Read lines from f until a blank line is encountered."""
593 while True:
594 line = f.readline()
595 if not line:
596 raise SpacegroupNotFoundError(
597 'invalid spacegroup `%s`, setting `%s` not found in data base'
598 % (spacegroup, setting))
599 if not line.strip():
600 break
603def _skip_to_nonblank(f, spacegroup, setting):
604 """Read lines from f until a nonblank line not starting with a
605 hash (#) is encountered and returns this and the next line."""
606 while True:
607 line1 = f.readline()
608 if not line1:
609 raise SpacegroupNotFoundError(
610 'invalid spacegroup %s, setting %i not found in data base' %
611 (spacegroup, setting))
612 line1.strip()
613 if line1 and not line1.startswith('#'):
614 line2 = f.readline()
615 break
616 return line1, line2
619def _read_datafile_entry(spg, no, symbol, setting, f):
620 """Read space group data from f to spg."""
622 floats = {'0.0': 0.0, '1.0': 1.0, '0': 0.0, '1': 1.0, '-1': -1.0}
623 for n, d in [(1, 2), (1, 3), (2, 3), (1, 4), (3, 4), (1, 6), (5, 6)]:
624 floats['{0}/{1}'.format(n, d)] = n / d
625 floats['-{0}/{1}'.format(n, d)] = -n / d
627 spg._no = no
628 spg._symbol = symbol.strip()
629 spg._setting = setting
630 spg._centrosymmetric = bool(int(f.readline().split()[1]))
631 # primitive vectors
632 f.readline()
633 spg._scaled_primitive_cell = np.array(
634 [[float(floats.get(s, s)) for s in f.readline().split()]
635 for i in range(3)],
636 dtype=float)
637 # primitive reciprocal vectors
638 f.readline()
639 spg._reciprocal_cell = np.array([[int(i) for i in f.readline().split()]
640 for i in range(3)],
641 dtype=int)
642 # subtranslations
643 spg._nsubtrans = int(f.readline().split()[0])
644 spg._subtrans = np.array(
645 [[float(floats.get(t, t)) for t in f.readline().split()]
646 for i in range(spg._nsubtrans)],
647 dtype=float)
648 # symmetry operations
649 nsym = int(f.readline().split()[0])
650 symop = np.array([[float(floats.get(s, s)) for s in f.readline().split()]
651 for i in range(nsym)],
652 dtype=float)
653 spg._nsymop = nsym
654 spg._rotations = np.array(symop[:, :9].reshape((nsym, 3, 3)), dtype=int)
655 spg._translations = symop[:, 9:]
658def _read_datafile(spg, spacegroup, setting, f):
659 if isinstance(spacegroup, int):
660 pass
661 elif isinstance(spacegroup, str):
662 spacegroup = ' '.join(spacegroup.strip().split())
663 compact_spacegroup = ''.join(spacegroup.split())
664 else:
665 raise SpacegroupValueError('`spacegroup` must be of type int or str')
666 while True:
667 line1, line2 = _skip_to_nonblank(f, spacegroup, setting)
668 _no, _symbol = line1.strip().split(None, 1)
669 _symbol = format_symbol(_symbol)
670 compact_symbol = ''.join(_symbol.split())
671 _setting = int(line2.strip().split()[1])
672 _no = int(_no)
673 if ((isinstance(spacegroup, int) and _no == spacegroup
674 and _setting == setting)
675 or (isinstance(spacegroup, str)
676 and compact_symbol == compact_spacegroup) and
677 (setting is None or _setting == setting)):
678 _read_datafile_entry(spg, _no, _symbol, _setting, f)
679 break
680 else:
681 _skip_to_blank(f, spacegroup, setting)
684def parse_sitesym_element(element):
685 """Parses one element from a single site symmetry in the form used
686 by the International Tables.
688 Examples:
690 >>> parse_sitesym_element("x")
691 ([(0, 1)], 0.0)
692 >>> parse_sitesym_element("-1/2-y")
693 ([(1, -1)], -0.5)
694 >>> parse_sitesym_element("z+0.25")
695 ([(2, 1)], 0.25)
696 >>> parse_sitesym_element("x-z+0.5")
697 ([(0, 1), (2, -1)], 0.5)
701 Parameters
702 ----------
704 element: str
705 Site symmetry like "x" or "-y+1/4" or "0.5+z".
708 Returns
709 -------
711 list[tuple[int, int]]
712 Rotation information in the form '(index, sign)' where index is
713 0 for "x", 1 for "y" and 2 for "z" and sign is '1' for a positive
714 entry and '-1' for a negative entry. E.g. "x" is '(0, 1)' and
715 "-z" is (2, -1).
717 float
718 Translation information in fractional space. E.g. "-1/4" is
719 '-0.25' and "1/2" is '0.5' and "0.75" is '0.75'.
722 """
723 element = element.lower()
724 is_positive = True
725 is_frac = False
726 sng_trans = None
727 fst_trans = []
728 snd_trans = []
729 rot = []
731 for char in element:
732 if char == "+":
733 is_positive = True
734 elif char == "-":
735 is_positive = False
736 elif char == "/":
737 is_frac = True
738 elif char in "xyz":
739 rot.append((ord(char)-ord("x"), 1 if is_positive else -1))
740 elif char.isdigit() or char == ".":
741 if sng_trans is None:
742 sng_trans = 1.0 if is_positive else -1.0
743 if is_frac:
744 snd_trans.append(char)
745 else:
746 fst_trans.append(char)
748 trans = 0.0 if not fst_trans else (sng_trans * float("".join(fst_trans)))
749 if is_frac:
750 trans /= float("".join(snd_trans))
752 return rot, trans
755def parse_sitesym_single(sym, out_rot, out_trans, sep=",", force_positive_translation=False):
756 """Parses a single site symmetry in the form used by International
757 Tables and overwrites 'out_rot' and 'out_trans' with data.
759 Parameters
760 ----------
762 sym: str
763 Site symmetry in the form used by International Tables (e.g. "x,y,z", "y-1/2,x,-z").
765 out_rot: np.array
766 A 3x3-integer array representing rotations (changes are made inplace).
768 out_rot: np.array
769 A 3-float array representing translations (changes are made inplace).
771 sep: str
772 String separator ("," in "x,y,z").
774 force_positive_translation: bool
775 Forces fractional translations to be between 0 and 1 (otherwise negative values might be accepted).
776 Defaults to 'False'.
779 Returns
780 -------
782 Nothing is returned: 'out_rot' and 'out_trans' are changed inplace.
785 """
786 out_rot[:] = 0.0
787 out_trans[:] = 0.0
789 for i, element in enumerate(sym.split(sep)):
790 e_rot_list, e_trans = parse_sitesym_element(element)
791 for rot_idx, rot_sgn in e_rot_list:
792 out_rot[i][rot_idx] = rot_sgn
793 out_trans[i] = (e_trans % 1.0) if force_positive_translation else e_trans
796def parse_sitesym(symlist, sep=',', force_positive_translation=False):
797 """Parses a sequence of site symmetries in the form used by
798 International Tables and returns corresponding rotation and
799 translation arrays.
801 Example:
803 >>> symlist = [
804 ... 'x,y,z',
805 ... '-y+1/2,x+1/2,z',
806 ... '-y,-x,-z',
807 ... 'x-1/4, y-1/4, -z'
808 ... ]
809 >>> rot, trans = parse_sitesym(symlist)
810 >>> rot
811 array([[[ 1, 0, 0],
812 [ 0, 1, 0],
813 [ 0, 0, 1]],
814 <BLANKLINE>
815 [[ 0, -1, 0],
816 [ 1, 0, 0],
817 [ 0, 0, 1]],
818 <BLANKLINE>
819 [[ 0, -1, 0],
820 [-1, 0, 0],
821 [ 0, 0, -1]],
822 <BLANKLINE>
823 [[ 1, 0, 0],
824 [ 0, 1, 0],
825 [ 0, 0, -1]]])
826 >>> trans
827 array([[ 0. , 0. , 0. ],
828 [ 0.5 , 0.5 , 0. ],
829 [ 0. , 0. , 0. ],
830 [-0.25, -0.25, 0. ]])
831 """
833 nsym = len(symlist)
834 rot = np.zeros((nsym, 3, 3), dtype='int')
835 trans = np.zeros((nsym, 3))
837 for i, sym in enumerate(symlist):
838 parse_sitesym_single(sym, rot[i], trans[i], sep=sep, force_positive_translation=force_positive_translation)
840 return rot, trans
843def spacegroup_from_data(no=None,
844 symbol=None,
845 setting=None,
846 centrosymmetric=None,
847 scaled_primitive_cell=None,
848 reciprocal_cell=None,
849 subtrans=None,
850 sitesym=None,
851 rotations=None,
852 translations=None,
853 datafile=None):
854 """Manually create a new space group instance. This might be
855 useful when reading crystal data with its own spacegroup
856 definitions."""
857 if no is not None and setting is not None:
858 spg = Spacegroup(no, setting, datafile)
859 elif symbol is not None:
860 spg = Spacegroup(symbol, None, datafile)
861 else:
862 raise SpacegroupValueError('either *no* and *setting* '
863 'or *symbol* must be given')
864 if not isinstance(sitesym, list):
865 raise TypeError('sitesym must be a list')
867 have_sym = False
868 if centrosymmetric is not None:
869 spg._centrosymmetric = bool(centrosymmetric)
870 if scaled_primitive_cell is not None:
871 spg._scaled_primitive_cell = np.array(scaled_primitive_cell)
872 if reciprocal_cell is not None:
873 spg._reciprocal_cell = np.array(reciprocal_cell)
874 if subtrans is not None:
875 spg._subtrans = np.atleast_2d(subtrans)
876 spg._nsubtrans = spg._subtrans.shape[0]
877 if sitesym is not None:
878 spg._rotations, spg._translations = parse_sitesym(sitesym)
879 have_sym = True
880 if rotations is not None:
881 spg._rotations = np.atleast_3d(rotations)
882 have_sym = True
883 if translations is not None:
884 spg._translations = np.atleast_2d(translations)
885 have_sym = True
886 if have_sym:
887 if spg._rotations.shape[0] != spg._translations.shape[0]:
888 raise SpacegroupValueError('inconsistent number of rotations and '
889 'translations')
890 spg._nsymop = spg._rotations.shape[0]
891 return spg
894def get_spacegroup(atoms, symprec=1e-5):
895 """Determine the spacegroup to which belongs the Atoms object.
897 This requires spglib: https://atztogo.github.io/spglib/ .
899 Parameters:
901 atoms: Atoms object
902 Types, positions and unit-cell.
903 symprec: float
904 Symmetry tolerance, i.e. distance tolerance in Cartesian
905 coordinates to find crystal symmetry.
907 The Spacegroup object is returned.
908 """
910 # Example:
911 # (We don't include the example in docstring to appease doctests
912 # when import fails)
913 # >>> from ase.build import bulk
914 # >>> atoms = bulk("Cu", "fcc", a=3.6, cubic=True)
915 # >>> sg = get_spacegroup(atoms)
916 # >>> sg
917 # Spacegroup(225, setting=1)
918 # >>> sg.no
919 # 225
921 import spglib
923 sg = spglib.get_spacegroup((atoms.get_cell(), atoms.get_scaled_positions(),
924 atoms.get_atomic_numbers()),
925 symprec=symprec)
926 if sg is None:
927 raise RuntimeError('Spacegroup not found')
928 sg_no = int(sg[sg.find('(') + 1:sg.find(')')])
929 return Spacegroup(sg_no)