Coverage for /builds/ase/ase/ase/spacegroup/spacegroup.py : 92.17%
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)
674 condition = (
675 (isinstance(spacegroup, int) and _no == spacegroup
676 and _setting == setting)
677 or (isinstance(spacegroup, str)
678 and compact_symbol == compact_spacegroup) and
679 (setting is None or _setting == setting))
681 if condition:
682 _read_datafile_entry(spg, _no, _symbol, _setting, f)
683 break
684 else:
685 _skip_to_blank(f, spacegroup, setting)
688def parse_sitesym_element(element):
689 """Parses one element from a single site symmetry in the form used
690 by the International Tables.
692 Examples:
694 >>> parse_sitesym_element("x")
695 ([(0, 1)], 0.0)
696 >>> parse_sitesym_element("-1/2-y")
697 ([(1, -1)], -0.5)
698 >>> parse_sitesym_element("z+0.25")
699 ([(2, 1)], 0.25)
700 >>> parse_sitesym_element("x-z+0.5")
701 ([(0, 1), (2, -1)], 0.5)
705 Parameters
706 ----------
708 element: str
709 Site symmetry like "x" or "-y+1/4" or "0.5+z".
712 Returns
713 -------
715 list[tuple[int, int]]
716 Rotation information in the form '(index, sign)' where index is
717 0 for "x", 1 for "y" and 2 for "z" and sign is '1' for a positive
718 entry and '-1' for a negative entry. E.g. "x" is '(0, 1)' and
719 "-z" is (2, -1).
721 float
722 Translation information in fractional space. E.g. "-1/4" is
723 '-0.25' and "1/2" is '0.5' and "0.75" is '0.75'.
726 """
727 element = element.lower()
728 is_positive = True
729 is_frac = False
730 sng_trans = None
731 fst_trans = []
732 snd_trans = []
733 rot = []
735 for char in element:
736 if char == "+":
737 is_positive = True
738 elif char == "-":
739 is_positive = False
740 elif char == "/":
741 is_frac = True
742 elif char in "xyz":
743 rot.append((ord(char) - ord("x"), 1 if is_positive else -1))
744 elif char.isdigit() or char == ".":
745 if sng_trans is None:
746 sng_trans = 1.0 if is_positive else -1.0
747 if is_frac:
748 snd_trans.append(char)
749 else:
750 fst_trans.append(char)
752 trans = 0.0 if not fst_trans else (sng_trans * float("".join(fst_trans)))
753 if is_frac:
754 trans /= float("".join(snd_trans))
756 return rot, trans
759def parse_sitesym_single(sym, out_rot, out_trans, sep=",",
760 force_positive_translation=False):
761 """Parses a single site symmetry in the form used by International
762 Tables and overwrites 'out_rot' and 'out_trans' with data.
764 Parameters
765 ----------
767 sym: str
768 Site symmetry in the form used by International Tables
769 (e.g. "x,y,z", "y-1/2,x,-z").
771 out_rot: np.array
772 A 3x3-integer array representing rotations (changes are made inplace).
774 out_rot: np.array
775 A 3-float array representing translations (changes are made inplace).
777 sep: str
778 String separator ("," in "x,y,z").
780 force_positive_translation: bool
781 Forces fractional translations to be between 0 and 1 (otherwise
782 negative values might be accepted). Defaults to 'False'.
785 Returns
786 -------
788 Nothing is returned: 'out_rot' and 'out_trans' are changed inplace.
791 """
792 out_rot[:] = 0.0
793 out_trans[:] = 0.0
795 for i, element in enumerate(sym.split(sep)):
796 e_rot_list, e_trans = parse_sitesym_element(element)
797 for rot_idx, rot_sgn in e_rot_list:
798 out_rot[i][rot_idx] = rot_sgn
799 out_trans[i] = \
800 (e_trans % 1.0) if force_positive_translation else e_trans
803def parse_sitesym(symlist, sep=',', force_positive_translation=False):
804 """Parses a sequence of site symmetries in the form used by
805 International Tables and returns corresponding rotation and
806 translation arrays.
808 Example:
810 >>> symlist = [
811 ... 'x,y,z',
812 ... '-y+1/2,x+1/2,z',
813 ... '-y,-x,-z',
814 ... 'x-1/4, y-1/4, -z'
815 ... ]
816 >>> rot, trans = parse_sitesym(symlist)
817 >>> rot
818 array([[[ 1, 0, 0],
819 [ 0, 1, 0],
820 [ 0, 0, 1]],
821 <BLANKLINE>
822 [[ 0, -1, 0],
823 [ 1, 0, 0],
824 [ 0, 0, 1]],
825 <BLANKLINE>
826 [[ 0, -1, 0],
827 [-1, 0, 0],
828 [ 0, 0, -1]],
829 <BLANKLINE>
830 [[ 1, 0, 0],
831 [ 0, 1, 0],
832 [ 0, 0, -1]]])
833 >>> trans
834 array([[ 0. , 0. , 0. ],
835 [ 0.5 , 0.5 , 0. ],
836 [ 0. , 0. , 0. ],
837 [-0.25, -0.25, 0. ]])
838 """
840 nsym = len(symlist)
841 rot = np.zeros((nsym, 3, 3), dtype='int')
842 trans = np.zeros((nsym, 3))
844 for i, sym in enumerate(symlist):
845 parse_sitesym_single(
846 sym, rot[i], trans[i], sep=sep,
847 force_positive_translation=force_positive_translation)
849 return rot, trans
852def spacegroup_from_data(no=None,
853 symbol=None,
854 setting=None,
855 centrosymmetric=None,
856 scaled_primitive_cell=None,
857 reciprocal_cell=None,
858 subtrans=None,
859 sitesym=None,
860 rotations=None,
861 translations=None,
862 datafile=None):
863 """Manually create a new space group instance. This might be
864 useful when reading crystal data with its own spacegroup
865 definitions."""
866 if no is not None and setting is not None:
867 spg = Spacegroup(no, setting, datafile)
868 elif symbol is not None:
869 spg = Spacegroup(symbol, None, datafile)
870 else:
871 raise SpacegroupValueError('either *no* and *setting* '
872 'or *symbol* must be given')
873 if not isinstance(sitesym, list):
874 raise TypeError('sitesym must be a list')
876 have_sym = False
877 if centrosymmetric is not None:
878 spg._centrosymmetric = bool(centrosymmetric)
879 if scaled_primitive_cell is not None:
880 spg._scaled_primitive_cell = np.array(scaled_primitive_cell)
881 if reciprocal_cell is not None:
882 spg._reciprocal_cell = np.array(reciprocal_cell)
883 if subtrans is not None:
884 spg._subtrans = np.atleast_2d(subtrans)
885 spg._nsubtrans = spg._subtrans.shape[0]
886 if sitesym is not None:
887 spg._rotations, spg._translations = parse_sitesym(sitesym)
888 have_sym = True
889 if rotations is not None:
890 spg._rotations = np.atleast_3d(rotations)
891 have_sym = True
892 if translations is not None:
893 spg._translations = np.atleast_2d(translations)
894 have_sym = True
895 if have_sym:
896 if spg._rotations.shape[0] != spg._translations.shape[0]:
897 raise SpacegroupValueError('inconsistent number of rotations and '
898 'translations')
899 spg._nsymop = spg._rotations.shape[0]
900 return spg
903def get_spacegroup(atoms, symprec=1e-5):
904 """Determine the spacegroup to which belongs the Atoms object.
906 This requires spglib: https://atztogo.github.io/spglib/ .
908 Parameters:
910 atoms: Atoms object
911 Types, positions and unit-cell.
912 symprec: float
913 Symmetry tolerance, i.e. distance tolerance in Cartesian
914 coordinates to find crystal symmetry.
916 The Spacegroup object is returned.
917 """
919 # Example:
920 # (We don't include the example in docstring to appease doctests
921 # when import fails)
922 # >>> from ase.build import bulk
923 # >>> atoms = bulk("Cu", "fcc", a=3.6, cubic=True)
924 # >>> sg = get_spacegroup(atoms)
925 # >>> sg
926 # Spacegroup(225, setting=1)
927 # >>> sg.no
928 # 225
930 import spglib
932 sg = spglib.get_spacegroup((atoms.get_cell(), atoms.get_scaled_positions(),
933 atoms.get_atomic_numbers()),
934 symprec=symprec)
935 if sg is None:
936 raise RuntimeError('Spacegroup not found')
937 sg_no = int(sg[sg.find('(') + 1:sg.find(')')])
938 return Spacegroup(sg_no)