Source code for ase.geometry.minkowski_reduction

import itertools
import numpy as np
from ase.utils import pbc2pbc
from ase.cell import Cell


TOL = 1E-12
MAX_IT = 100000    # in practice this is not exceeded


class CycleChecker:

    def __init__(self, d):
        assert d in [2, 3]

        # worst case is the hexagonal cell in 2D and the fcc cell in 3D
        n = {2: 6, 3: 12}[d]

        # max cycle length is total number of primtive cell descriptions
        max_cycle_length = np.prod([n - i for i in range(d)]) * np.prod(d)
        self.visited = np.zeros((max_cycle_length, 3 * d), dtype=int)

    def add_site(self, H):
        # flatten array for simplicity
        H = H.ravel()

        # check if site exists
        found = (self.visited == H).all(axis=1).any()

        # shift all visited sites down and place current site at the top
        self.visited = np.roll(self.visited, 1, axis=0)
        self.visited[0] = H
        return found


def reduction_gauss(B, hu, hv):
    """Calculate a Gauss-reduced lattice basis (2D reduction)."""
    cycle_checker = CycleChecker(d=2)
    u = hu @ B
    v = hv @ B

    for it in range(MAX_IT):
        x = int(round(np.dot(u, v) / np.dot(u, u)))
        hu, hv = hv - x * hu, hu
        u = hu @ B
        v = hv @ B
        site = np.array([hu, hv])
        if np.dot(u, u) >= np.dot(v, v) or cycle_checker.add_site(site):
            return hv, hu

    raise RuntimeError(f"Gaussian basis not found after {MAX_IT} iterations")


def relevant_vectors_2D(u, v):
    cs = np.array([e for e in itertools.product([-1, 0, 1], repeat=2)])
    vs = cs @ [u, v]
    indices = np.argsort(np.linalg.norm(vs, axis=1))[:7]
    return vs[indices], cs[indices]


def closest_vector(t0, u, v):
    t = t0
    a = np.zeros(2, dtype=int)
    rs, cs = relevant_vectors_2D(u, v)

    dprev = float("inf")
    for it in range(MAX_IT):
        ds = np.linalg.norm(rs + t, axis=1)
        index = np.argmin(ds)
        if index == 0 or ds[index] >= dprev:
            return a

        dprev = ds[index]
        r = rs[index]
        kopt = int(round(-np.dot(t, r) / np.dot(r, r)))
        a += kopt * cs[index]
        t = t0 + a[0] * u + a[1] * v

    raise RuntimeError(f"Closest vector not found after {MAX_IT} iterations")


def reduction_full(B):
    """Calculate a Minkowski-reduced lattice basis (3D reduction)."""
    cycle_checker = CycleChecker(d=3)
    H = np.eye(3, dtype=int)
    norms = np.linalg.norm(B, axis=1)

    for it in range(MAX_IT):
        # Sort vectors by norm
        H = H[np.argsort(norms, kind='merge')]

        # Gauss-reduce smallest two vectors
        hw = H[2]
        hu, hv = reduction_gauss(B, H[0], H[1])
        H = np.array([hu, hv, hw])
        R = H @ B

        # Orthogonalize vectors using Gram-Schmidt
        u, v, _ = R
        X = u / np.linalg.norm(u)
        Y = v - X * np.dot(v, X)
        Y /= np.linalg.norm(Y)

        # Find closest vector to last element of R
        pu, pv, pw = R @ np.array([X, Y]).T
        nb = closest_vector(pw, pu, pv)

        # Update basis
        H[2] = [nb[0], nb[1], 1] @ H
        R = H @ B

        norms = np.linalg.norm(R, axis=1)
        if norms[2] >= norms[1] or cycle_checker.add_site(H):
            return R, H

    raise RuntimeError(f"Reduced basis not found after {MAX_IT} iterations")


[docs]def is_minkowski_reduced(cell, pbc=True): """Tests if a cell is Minkowski-reduced. Parameters: cell: array The lattice basis to test (in row-vector format). pbc: array, optional The periodic boundary conditions of the cell (Default `True`). If `pbc` is provided, only periodic cell vectors are tested. Returns: is_reduced: bool True if cell is Minkowski-reduced, False otherwise. """ """These conditions are due to Minkowski, but a nice description in English can be found in the thesis of Carine Jaber: "Algorithmic approaches to Siegel's fundamental domain", https://www.theses.fr/2017UBFCK006.pdf This is also good background reading for Minkowski reduction. 0D and 1D cells are trivially reduced. For 2D cells, the conditions which an already-reduced basis fulfil are: |b1| ≤ |b2| |b2| ≤ |b1 - b2| |b2| ≤ |b1 + b2| For 3D cells, the conditions which an already-reduced basis fulfil are: |b1| ≤ |b2| ≤ |b3| |b1 + b2| ≥ |b2| |b1 + b3| ≥ |b3| |b2 + b3| ≥ |b3| |b1 - b2| ≥ |b2| |b1 - b3| ≥ |b3| |b2 - b3| ≥ |b3| |b1 + b2 + b3| ≥ |b3| |b1 - b2 + b3| ≥ |b3| |b1 + b2 - b3| ≥ |b3| |b1 - b2 - b3| ≥ |b3| """ pbc = pbc2pbc(pbc) dim = pbc.sum() if dim <= 1: return True if dim == 2: # reorder cell vectors to [shortest, longest, aperiodic] cell = cell.copy() cell[np.argmin(pbc)] = 0 norms = np.linalg.norm(cell, axis=1) cell = cell[np.argsort(norms)[[1, 2, 0]]] A = [[0, 1, 0], [1, -1, 0], [1, 1, 0]] lhs = np.linalg.norm(A @ cell, axis=1) norms = np.linalg.norm(cell, axis=1) rhs = norms[[0, 1, 1]] else: A = [[0, 1, 0], [0, 0, 1], [1, 1, 0], [1, 0, 1], [0, 1, 1], [1, -1, 0], [1, 0, -1], [0, 1, -1], [1, 1, 1], [1, -1, 1], [1, 1, -1], [1, -1, -1]] lhs = np.linalg.norm(A @ cell, axis=1) norms = np.linalg.norm(cell, axis=1) rhs = norms[[0, 1, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2]] return (lhs >= rhs - TOL).all()
[docs]def minkowski_reduce(cell, pbc=True): """Calculate a Minkowski-reduced lattice basis. The reduced basis has the shortest possible vector lengths and has norm(a) <= norm(b) <= norm(c). Implements the method described in: Low-dimensional Lattice Basis Reduction Revisited Nguyen, Phong Q. and Stehlé, Damien, ACM Trans. Algorithms 5(4) 46:1--46:48, 2009 :doi:`10.1145/1597036.1597050` Parameters: cell: array The lattice basis to reduce (in row-vector format). pbc: array, optional The periodic boundary conditions of the cell (Default `True`). If `pbc` is provided, only periodic cell vectors are reduced. Returns: rcell: array The reduced lattice basis. op: array The unimodular matrix transformation (rcell = op @ cell). """ cell = Cell(cell) pbc = pbc2pbc(pbc) dim = pbc.sum() op = np.eye(3, dtype=int) if is_minkowski_reduced(cell, pbc): return cell, op if dim == 2: # permute cell so that first two vectors are the periodic ones perm = np.argsort(pbc, kind='merge')[::-1] # stable sort pcell = cell[perm][:, perm] # perform gauss reduction norms = np.linalg.norm(pcell, axis=1) norms[2] = float("inf") indices = np.argsort(norms) op = op[indices] hu, hv = reduction_gauss(pcell, op[0], op[1]) op[0] = hu op[1] = hv # undo above permutation invperm = np.argsort(perm) op = op[invperm][:, invperm] # maintain cell handedness index = np.argmin(pbc) normal = np.cross(cell[index - 2], cell[index - 1]) normal /= np.linalg.norm(normal) _cell = cell.copy() _cell[index] = normal _rcell = op @ cell _rcell[index] = normal if _cell.handedness != Cell(_rcell).handedness: op[index - 1] *= -1 elif dim == 3: _, op = reduction_full(cell) # maintain cell handedness if cell.handedness != Cell(op @ cell).handedness: op = -op norms1 = np.sort(np.linalg.norm(cell, axis=1)) norms2 = np.sort(np.linalg.norm(op @ cell, axis=1)) if (norms2 > norms1 + TOL).any(): raise RuntimeError("Minkowski reduction failed") return op @ cell, op