niggli.py

"""
This module finds the Niggli cell representation of a given unit cell.
The Niggli cell is a uniquely-defined maximally-reduced unit cell.
Relevant citations are

Niggli, P. "Krystallographische und strukturtheoretische Grundbegriffe.
Handbuch der Experimentalphysik", 1928, Vol. 7, Part 1, 108-176.

Krivy, I. and Gruber, B., "A Unified Algorithm for Determining the 
Reduced (Niggli) Cell", Acta Cryst. 1976, A32, 297-298.

Grosse-Kunstleve, R.W.; Sauter, N. K.; and Adams, P. D. "Numerically
stable algorithms for the computation of reduced unit cells", Acta Cryst.
2004, A60, 1-6.
"""

import numpy as np
from ase.lattice.spacegroup.cell import cellpar_to_cell


class _gtensor(object):
    """
    The G tensor as defined in Grosse-Kunstleve.
    """
    def __init__(self, cell, epsilon):
        self.cell = cell
        self.epsilon = epsilon

        self.a = np.dot(cell[0], cell[0])
        self.b = np.dot(cell[1], cell[1])
        self.c = np.dot(cell[2], cell[2])

        self.x = 2 * np.dot(cell[1], cell[2])
        self.y = 2 * np.dot(cell[0], cell[2])
        self.z = 2 * np.dot(cell[0], cell[1])

        self._G = np.array([
            [self.a, self.z/2., self.y/2.],
            [self.z/2., self.b, self.x/2.],
            [self.y/2., self.x/2., self.c],
            ])

        self._lmn()

    def update(self, C):
        """
        Procedure A0 as defined in Krivy.
        """
        self._G = np.dot(C.T, np.dot(self._G, C))

        self.a = self._G[0][0]
        self.b = self._G[1][1]
        self.c = self._G[2][2]
    
        self.x = 2 * self._G[1][2]
        self.y = 2 * self._G[0][2]
        self.z = 2 * self._G[0][1]

        self._lmn()

    def _lmn(self):
        self.l = 0
        self.m = 0
        self.n = 0

        if self.x < -self.epsilon:
            self.l = -1
        elif self.x > self.epsilon:
            self.l = 1
        if self.y < -self.epsilon:
            self.m = -1
        elif self.y > self.epsilon:
            self.m = 1
        if self.z < -self.epsilon:
            self.n = -1
        elif self.z > self.epsilon:
            self.n = 1

    def get_newcell(self):
        a = np.sqrt(self.a)
        b = np.sqrt(self.b)
        c = np.sqrt(self.c)

        alpha = np.arccos(self.x / (2 * b * c)) * 180. / np.pi
        beta = np.arccos(self.y / (2 * a * c)) * 180. / np.pi
        gamma = np.arccos(self.z / (2 * a * b)) * 180. / np.pi

        return cellpar_to_cell([a, b, c, alpha, beta, gamma])


def niggli(cell):
    """ 
    This procedure takes a unit cell and returns its fully reduced Niggli
    representation.
    """
    V = np.dot(cell[0], np.cross(cell[1], cell[2]))
    G = _gtensor(cell, 1e-5 * V**(1./3.))

    # Once A2 and A5-A8 all evaluate to False, the unit cell will have
    # been fully reduced.
    for count in xrange(10000):
        if (G.a > G.b + G.epsilon or
                (not np.abs(G.a - G.b) > G.epsilon
                    and np.abs(G.x) > np.abs(G.y) + G.epsilon)):
            # Procedure A1
            G.update(np.array([
                [0, -1, 0],
                [-1, 0, 0],
                [0, 0, -1],
                ]))

        if (G.b > G.c + G.epsilon or
                (not np.abs(G.b - G.c) > G.epsilon
                    and np.abs(G.y) > np.abs(G.z) + G.epsilon)):
            # Procedure A2
            G.update(np.array([
                [-1, 0, 0],
                [0, 0, -1],
                [0, -1, 0],
                ]))
            continue

        if G.l * G.m * G.n == 1:
            # Procedure A3
            i = -1 if G.l == -1 else 1
            j = -1 if G.m == -1 else 1
            k = -1 if G.n == -1 else 1
            G.update(np.array([
                [i, 0, 0],
                [0, j, 0],
                [0, 0, k],
                ]))
        else:
            # Procedure A4
            i = -1 if G.l == 1 else 1
            j = -1 if G.m == 1 else 1
            k = -1 if G.n == 1 else 1

            if (i * j * k == -1):
                if G.l == 0:
                    i = -1
                if G.m == 0:
                    j = -1
                if G.n == 0:
                    k = -1
            G.update(np.array([
                [i, 0, 0],
                [0, j, 0],
                [0, 0, k],
                ]))

        if (np.abs(G.x) > G.b + G.epsilon
                or (not np.abs(G.b - G.x) > G.epsilon
                    and 2 * G.y < G.z - G.epsilon)
                or (not np.abs(G.b + G.x) > G.epsilon
                    and G.z < -G.epsilon)):
            # Procedure A5
            G.update(np.array([
                [1, 0, 0],
                [0, 1, -np.sign(G.x)],
                [0, 0, 1],
                ]))
        elif (np.abs(G.y) > G.a + G.epsilon
                or (not np.abs(G.a - G.y) > G.epsilon
                    and 2 * G.x < G.z - G.epsilon)
                or (not np.abs(G.a + G.y) > G.epsilon
                    and G.z < -G.epsilon)):
            # Procedure A6
            G.update(np.array([
                [1, 0, -np.sign(G.y)],
                [0, 1, 0],
                [0, 0, 1],
                ]))
        elif (np.abs(G.z) > G.a + G.epsilon
                or (not np.abs(G.a - G.z) > G.epsilon
                    and 2 * G.x < G.y - G.epsilon)
                or (not np.abs(G.a + G.z) > G.epsilon
                    and G.y < -G.epsilon)):
            # Procedure A7
            G.update(np.array([
                [1, -np.sign(G.z), 0],
                [0, 1, 0],
                [0, 0, 1],
                ]))
        elif (G.x + G.y + G.z + G.a + G.b < -G.epsilon
                or (not np.abs(G.x + G.y + G.z + G.a + G.b) > G.epsilon
                    and 2*(G.a + G.y) + G.z > G.epsilon)):
            # Procedure A8
            G.update(np.array([
                [1, 0, 1],
                [0, 1, 1],
                [0, 0, 1],
                ]))
        else:
            break
    else:
        raise RuntimeError('Niggli did not converge \
                in {n} iterations!'.format(n=count))
    return G.get_newcell()