Chimie Théorique » QMX

# Copyright (C) 2010, Jesper Friis

# (see accompanying license files for details).

"""Definition of the Spacegroup class.

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

"""

import os

import warnings

import numpy as np

__all__ = ['Spacegroup']

class SpacegroupError(Exception):

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

    pass

class SpacegroupNotFoundError(SpacegroupError):

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

    pass

class SpacegroupValueError(SpacegroupError):

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

    pass

class Spacegroup(object):

    """A space group class. 

    The instances of Spacegroup describes the symmetry operations for

    the given space group. 

    Example:

    >>> from spacegroup import Spacegroup

    >>> 

    >>> sg = Spacegroup(225)

    >>> print 'Space group', sg.no, sg.symbol

    Space group 225 F m -3 m

    >>> sg.scaled_primitive_cell

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

           [ 0.5,  0. ,  0.5],

           [ 0.5,  0.5,  0. ]])

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

    >>> sites

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

           [ 0. ,  0.5,  0.5],

           [ 0.5,  0. ,  0.5],

           [ 0.5,  0.5,  0. ]])

    """

    no = property(

        lambda self: self._no,

        doc='Space group number in International Tables of Crystallography.')

    symbol = property(

        lambda self: self._symbol, 

        doc='Hermann-Mauguin (or international) symbol for the space group.')

    setting = property(

        lambda self: self._setting, 

        doc='Space group setting. Either one or two.')

    lattice = property(

        lambda self: self._symbol[0], 

        doc="""Lattice type: 

            P      primitive

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

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

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

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

            """)

    centrosymmetric = property(

        lambda self: self._centrosymmetric, 

        doc='Whether a center of symmetry exists.')

    scaled_primitive_cell = property(

        lambda self: self._scaled_primitive_cell, 

        doc='Primitive cell in scaled coordinates as a matrix with the '

        'primitive vectors along the rows.')

    reciprocal_cell = property(

        lambda self: self._reciprocal_cell, 

        doc='Tree Miller indices that span all kinematically non-forbidden '

        'reflections as a matrix with the Miller indices along the rows.')

    nsubtrans = property(

        lambda self: len(self._subtrans), 

        doc='Number of cell-subtranslation vectors.')

    def _get_nsymop(self):

        if self.centrosymmetric:

            return 2 * len(self._rotations) * len(self._subtrans)

        else:

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

    nsymop = property(_get_nsymop, doc='Total number of symmetry operations.')

    subtrans = property(

        lambda self: self._subtrans, 

        doc='Translations vectors belonging to cell-sub-translations.')

    rotations = property(

        lambda self: self._rotations, 

        doc='Symmetry rotation matrices. The invertions are not included '

        'for centrosymmetrical crystals.')

    translations = property(

        lambda self: self._translations, 

        doc='Symmetry translations. The invertions are not included '

        'for centrosymmetrical crystals.')

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

        """Returns a new Spacegroup instance. 

        Parameters:

        spacegroup : int | string | Spacegroup instance

            The space group number in International Tables of

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

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

        setting : 1 | 2

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

            determines Which of these should be used.

        datafile : None | string

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

            will be used.

        """

        if isinstance(spacegroup, Spacegroup):

            for k, v in spacegroup.__dict__.iteritems():

                setattr(self, k, v)

            return 

        if not datafile:

            datafile = get_datafile()

        f = open(datafile, 'r')

        try:

            _read_datafile(self, spacegroup, setting, f)

        finally:

            f.close()

    def __repr__(self):

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

    def __str__(self):

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

        the same format as found the database."""

        retval = []

        # no, symbol

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

        # setting

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

        # centrosymmetric

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

        # primitive vectors

        retval.append('  primitive vectors\n')

        for i in range(3):

            retval.append('   ')

            for j in range(3):

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

            retval.append('\n')

        # primitive reciprocal vectors

        retval.append('  reciprocal vectors\n')

        for i in range(3):

            retval.append('   ')

            for j in range(3):

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

            retval.append('\n')

        # sublattice

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

        for i in range(self.nsubtrans):

            retval.append('   ')

            for j in range(3):

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

            retval.append('\n')

        # symmetry operations

        nrot = len(self.rotations)

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

        for i in range(nrot):

            retval.append(' ')

            for j in range(3):

                retval.append(' ')

                for k in range(3):

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

                retval.append('  ')

            for j in range(3):

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

            retval.append('\n')

        retval.append('\n')

        return ''.join(retval)

    def __eq__(self, other):

        """Chech whether *self* and *other* refer to the same

        spacegroup number and setting."""

        if isinstance(other, int):

            return self.no == other

        elif isinstance(other, str):

            return self.symbol == other

        else:

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

    def __index__(self):

        return self.no

    def get_symop(self):

        """Returns all symmetry operations (including inversions and

        subtranslations) as a sequence of (rotation, translation)

        tuples."""

        symop = []

        parities = [1]

        if self.centrosymmetric:

            parities.append(-1)

        for parity in parities:

            for subtrans in self.subtrans:

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

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

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

        return symop

    def get_rotations(self):

        """Return all rotations, including inversions for

        centrosymmetric crystals."""

        if self.centrosymmetric:

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

        else:

            return self.rotations

    def equivalent_reflections(self, hkl):

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

        in hkl.

        Example:

        >>> from spacegroup import Spacegroup

        >>> sg = Spacegroup(225)  # fcc

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

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

               [ 0, -2,  0],

               [-2,  0,  0],

               [ 2,  0,  0],

               [ 0,  2,  0],

               [ 0,  0,  2]])

        """

        hkl = np.array(hkl, dtype='int')

        rot = self.get_rotations()

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

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

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

        ind = np.lexsort(refl.T)

        refl = refl[ind]

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

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

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

    def equivalent_sites(self, scaled_positions, ondublicates='error', 

                         symprec=1e-3):

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

        Parameters:

        scaled_positions: list | array

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

        ondublicates : 'keep' | 'replace' | 'warn' | 'error'

            Action if `scaled_positions` contain symmetry-equivalent

            positions:

            'keep'

               ignore additional symmetry-equivalent positions

            'replace'

                replace

            'warn'

                like 'keep', but issue an UserWarning

            'error'

                raises a SpacegroupValueError

        symprec: float

            Minimum "distance" betweed two sites in scaled coordinates

            before they are counted as the same site.

        Returns:

        sites: array

            A NumPy array of equivalent sites.

        kinds: list

            A list of integer indices specifying which input site is 

            equivalent to the corresponding returned site.

        """

        kinds = []

        sites = []

        symprec2 = symprec**2

        for kind, pos in enumerate(scaled_positions):

            for rot, trans in self.get_symop():

                site = np.mod(np.dot(rot, pos) + trans, 1.)

                if not sites:

                    sites.append(site)

                    kinds.append(kind)

                    continue

                t = site - sites

                mask = np.sum(t*t, 1) < symprec2

                if np.any(mask):

                    ind = np.argwhere(mask)[0][0]

                    if kinds[ind] == kind:

                        pass

                    elif ondublicates == 'keep':

                        pass

                    elif ondublicates == 'replace':

                        kinds[ind] = kind

                    elif ondublicates == 'warn':

                        warnings.warn('scaled_positions %d and %d '

                                      'are equivalent'%(kinds[ind], kind))

                    elif ondublicates == 'error':

                        raise SpacegroupValueError(

                            'scaled_positions %d and %d are equivalent'%(

                                kinds[ind], kind))

                    else:

                        raise SpacegroupValueError(

                            'Argument "ondublicates" must be one of: '

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

                else:

                    sites.append(site)

                    kinds.append(kind)

        return np.array(sites), kinds

def get_datafile():

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

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

#-----------------------------------------------------------------

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

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

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

# instance is created.

#-----------------------------------------------------------------

def _skip_to_blank(f, spacegroup, setting):

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

    while True:

        line = f.readline()

        if not line:

            raise SpacegroupNotFoundError(

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

                ( spacegroup, setting ) )

        if not line.strip():

            break

def _skip_to_nonblank(f, spacegroup, setting):

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

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

    while True:

        line1 = f.readline()

        if not line1:

            raise SpacegroupNotFoundError(

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

                ( spacegroup, setting ) )

        line1.strip()

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

            line2 = f.readline()

            break

    return line1, line2

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

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

    spg._no = no

    spg._symbol = symbol.strip()

    spg._setting = setting

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

    # primitive vectors

    f.readline()

    spg._scaled_primitive_cell = np.array([map(float, f.readline().split()) 

                                       for i in range(3)], 

                                      dtype=np.float)

    # primitive reciprocal vectors

    f.readline()

    spg._reciprocal_cell = np.array([map(int, f.readline().split()) 

                                        for i in range(3)],

                                       dtype=np.int)

    # subtranslations

    spg._nsubtrans = int(f.readline().split()[0])

    spg._subtrans = np.array([map(float, f.readline().split()) 

                              for i in range(spg._nsubtrans)],

                             dtype=np.float)

    # symmetry operations

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

    symop = np.array([map(float, f.readline().split()) for i in range(nsym)],

                     dtype=np.float)

    spg._nsymop = nsym

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

    spg._translations = symop[:,9:]

def _read_datafile(spg, spacegroup, setting, f):

    if isinstance(spacegroup, int):

        pass

    elif isinstance(spacegroup, basestring):

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

    else:

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

    while True:

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

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

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

        _no = int(_no)

        if ((isinstance(spacegroup, int) and _no == spacegroup) or

            (isinstance(spacegroup, basestring) and 

             _symbol == spacegroup)) and _setting == setting:

            _read_datafile_entry(spg, _no, _symbol, _setting, f)

            break

        else:

            _skip_to_blank(f, spacegroup, setting)

#-----------------------------------------------------------------

# Self test

if __name__ == '__main__':

    # Import spacegroup in order to ensure that __file__ is defined

    # such that the data base can be found.

    import spacegroup

    import doctest

    print 'doctest: ', doctest.testmod()