Commit 691d8c85 authored by Dániel Varjas's avatar Dániel Varjas Committed by Joseph Weston

improve lll and voronoi algorithms in kwant.linalg

Improve the logic of 'cvp', now it provably works in all cases.
Add 'reduced' option to 'linalg.lll.voronoi', this eliminates
superfluous Voronoi vectors that do not have a corresponding face
of the Voronoi cell, or this face would be smaller than 'rtol'.

Simplify 'lattice.Polyatomic.neighbors' by providing the option
to group equidistant vectors 'group_by_length' in 'linalg.lll.cvp'.

See !165 for details.

Closes #163.
parent b416e486
Pipeline #12574 passed with stages
in 16 minutes and 24 seconds
......@@ -286,7 +286,8 @@ class Polyatomic:
n : integer
Order of the hoppings to return.
Order of the hoppings to return. Note that the zeroth neighbor is
the site itself or any other sites with the same position.
eps : float
Tolerance relative to the length of the shortest lattice vector for
when to consider lengths to be approximately equal.
......@@ -309,6 +310,7 @@ class Polyatomic:
sls = self.sublattices
shortest_hopping = sls[0].n_closest(
sls[0].pos(([0] * sls[0].lattice_dim)), 2)[-1]
rtol = eps
eps *= np.linalg.norm(self.vec(shortest_hopping))
nvec = len(self._prim_vecs)
sublat_pairs = [(i, j) for (i, j) in product(sls, sls)
......@@ -323,36 +325,25 @@ class Polyatomic:
return True
# Find the correct number of neighbors to calculate on each lattice.
cutoff = n + 2
while True:
max_dist = []
# Find the `n` closest neighbors (with multiplicity) for each
# pair of lattices, this surely includes the `n` closest neighbors overall.
sites = []
for i, j in sublat_pairs:
origin = Site(j, ta.zeros(nvec)).pos
tags = i.n_closest(origin, n=cutoff**nvec)
tags = i.n_closest(origin, n=n+1, group_by_length=True, rtol=rtol)
ij_dist = [np.linalg.norm(Site(i, tag).pos - origin)
for tag in tags]
sites.append((tags, (j, i), ij_dist))
max_dist = [i[2][-1] for i in sites]
distances = np.r_[tuple((i[2] for i in sites))]
distances = np.sort(distances)
group_boundaries = np.argwhere(np.diff(distances) > eps)
if len(group_boundaries) < n:
cutoff += 1
group_boundaries = np.where(np.diff(distances) > eps)[0]
# Find distance in `n`-th group.
if len(group_boundaries) == n:
n_dist = distances[-1]
n_dist = distances[group_boundaries[n]]
except IndexError:
cutoff += 1
if np.all(max_dist > n_dist):
cutoff += 1
# We now have all the required sites, we need to find n-th.
# We now have all the required sites, select the ones in `n`-th group.
result = []
for group in sites:
tags, distance = group[0], group[2]
......@@ -485,7 +476,7 @@ class Monatomic(builder.SiteFamily, Polyatomic):
raise ValueError("Dimensionality mismatch.")
return tag
def n_closest(self, pos, n=1):
def n_closest(self, pos, n=1, group_by_length=False, rtol=1e-9):
"""Find n sites closest to position `pos`.
......@@ -494,7 +485,8 @@ class Monatomic(builder.SiteFamily, Polyatomic):
An array with sites coordinates.
# TODO (Anton): transform to tinyarrays, once ta indexing is better.
return - self.offset, self.reduced_vecs, n),
return - self.offset, self.reduced_vecs,
n=n, group_by_length=group_by_length, rtol=rtol),
def closest(self, pos):
......@@ -9,6 +9,7 @@
__all__ = ['lll', 'cvp', 'voronoi']
import numpy as np
from scipy import linalg as la
from itertools import product
......@@ -109,12 +110,12 @@ def lll(basis, c=1.34):
return vecs, np.array(np.round(coefs), int)
def cvp(vec, basis, n=1):
def cvp(vec, basis, n=1, group_by_length=False, rtol=1e-09):
Solve the n-closest vector problem for a vector, given a basis.
This algorithm performs poorly in general, so it should be supplied
with LLL-reduced bases.
with LLL-reduced basis.
......@@ -123,12 +124,21 @@ def cvp(vec, basis, n=1):
basis : 2d array-like of floats
Sequence of basis vectors
n : int
Number of lattice vectors closest to the point that need to be found.
Number of lattice vectors closest to the point that need to be found. If
`group_by_length=True`, all vectors with the `n` shortest distinct distances
are found.
group_by_length : bool
Whether to count points with distances that are within `rtol` as one
towards `n`. Useful for finding `n`-th nearest neighbors of high symmetry points.
rtol : float
If `group_by_length=True`, relative tolerance when deciding whether
distances are equal.
coords : numpy array
An array with the coefficients of the lattice vectors closest to the
An array with the coefficients of the `n` closest lattice vectors, or all
the lattice vectors with the `n` shortest distances from the
requested point.
......@@ -142,32 +152,84 @@ def cvp(vec, basis, n=1):
if basis.ndim != 2:
raise ValueError('`basis` must be a 2d array-like object.')
vec = np.asarray(vec)
# TODO: change to rcond=None once we depend on numpy >= 1.14.
center_coords = np.linalg.lstsq(basis.T, vec, rcond=-1)[0]
center_coords = np.array(np.round(center_coords), int)
# Cutoff radius for n-th nearest neighbor.
rad = 1
nth_dist = np.inf
# Project the coordinates on the space spanned by `basis`, in
# units of `basis` vectors.
vec_bas = la.lstsq(basis.T, vec)[0]
# Round the coordinates, this finds the lattice point which is
# the center of the parallelepiped that contains `vec`.
# `vec` is surely in the parallelepiped with edges `basis`
# centered at `center_coords`.
center_coords = np.array(np.round(vec_bas), int)
# Make `vec` the projection of `vec` onto the space spanned by `basis`.
vec = vec_bas @ basis
bbt = basis @ basis.T
# The volume of a parallelepiped spanned by `basis` is `sqrt(det(bbt))`.
# (We use this instead of `det(basis)` because basis does not necessarily
# span the full space, hence not always a square matrix.)
# We calculate the radius `rad` of the largest sphere that can be inscribed
# in the parallelepiped spanned by `basis`. The area of each face spanned
# by one less vector in `basis` is given by the above formula applied to `bbt`
# with one row and column discarded (first minor). We utilizie the relation
# between entries of the inverse matrix and first minors. The perpendicular
# size of the parallelepiped is given by the volume of the parallelepiped
# divided by the area of the face. We choose the smallest perpendicular
# size as the diameter of the largest sphere inside.
rad = 0.5 / np.sqrt(np.max(np.diag(la.inv(bbt))))
l = 1
while True:
r = np.round(rad * np.linalg.cond(basis)) + 1
points = np.mgrid[tuple(slice(i - r, i + r) for i in center_coords)]
# Make all the lattice points in and on the edges of a parallelepiped
# of `2*l` size around `center_coords`.
points = np.mgrid[tuple(slice(i - l, i + l + 1) for i in center_coords)]
points = points.reshape(basis.shape[0], -1).T
# If there are less than `n` points, make more.
if len(points) < n:
rad += 1
l += 1
point_coords =, basis)
point_coords -= vec.T
distances = np.sqrt(np.sum(point_coords**2, 1))
point_coords = points @ basis
point_coords = point_coords - vec.T
distances = la.norm(point_coords, axis = 1)
order = np.argsort(distances)
new_nth_dist = distances[order[n - 1]]
if new_nth_dist < nth_dist:
nth_dist = new_nth_dist
rad += 1
distances = distances[order]
points = points[order]
if group_by_length:
# Group points that are equidistant within `rtol * rad`.
group_boundaries = np.where(np.diff(distances) > rtol * rad)[0]
if len(group_boundaries) + 1 < n:
# Make more points if there are less than `n` groups.
l += 1
elif len(group_boundaries) == n - 1:
# If there are exactly `n` groups, choose largest distance.
dist_n = distances[-1]
points_keep = points
# If there are more than `n` groups, choose the largest
# distance in the `n`-th group.
dist_n = distances[group_boundaries[n-1]]
# Only return points that are in the first `n` groups.
points_keep = points[:group_boundaries[n-1] + 1]
rtol = 0
# Choose the `n`-th distance.
dist_n = distances[n-1]
# Only return the first `n` points.
points_keep = points[:n]
# We covered all points within radius `(2*l-1) * rad` from the
# parallelepiped spanned by `basis` centered at `center_coords`.
# Because `vec` is guaranteed to be in this parallelepiped, if
# the current `dist_n` is smaller than `(2*l-1) * rad`, we surely
# found all the smaller vectors.
if dist_n < (2*l - 1 - rtol) * rad:
return np.array(points_keep, int)
return np.array(points[order[:n]], int)
# Otherwise there may be smaller vectors we haven't found,
# increase `l`.
l += 1
def voronoi(basis):
def voronoi(basis, reduced=False, rtol=1e-09):
Return an array of lattice vectors forming its voronoi cell.
......@@ -175,27 +237,53 @@ def voronoi(basis):
basis : 2d array-like of floats
Basis vectors for which the Voronoi neighbors have to be found.
reduced : bool
If False, exactly `2 (2^D - 1)` vectors are returned (where `D`
is the number of lattice vectors), these are not always the minimal
set of Voronoi vectors.
If True, only the minimal set of Voronoi vectors are returned.
rtol : float
Tolerance when deciding whether a vector is in the minimal set,
vectors associated with small faces compared to the size of the
unit cell are discarded.
voronoi_neighbors : numpy array of ints
All the lattice vectors that may potentially neighbor the origin.
This algorithm does not calculate the minimal Voronoi cell of the lattice
and can be optimized. Its main aim is flood-fill, however, and better
safe than sorry.
All the lattice vectors that (potentially) neighbor the origin.
List of length `2n` where the second half of the list contains
is `-1` times the vectors in the first half.
basis = np.asarray(basis)
if basis.ndim != 2:
raise ValueError('`basis` must be a 2d array-like object.')
# Find halved lattice points, every face of the VC contains half
# of the lattice vector which is the normal of the face.
# These points are all potentially on a face a VC,
# but not necessarily the VC centered at the origin.
displacements = list(product(*(len(basis) * [[0, .5]])))[1:]
vertices = np.array([cvp(, basis), basis)[0] for vec in
# Find the nearest lattice point, this is the lattice point whose
# VC this face belongs to.
vertices = np.array([cvp(vec @ basis, basis)[0] for vec in
# The lattice vector for a face is exactly twice the vector to the
# closest lattice point from the halved lattice point on the face.
vertices = np.array(np.round((vertices - displacements) * 2), int)
if reduced:
# Discard vertices that are not associated with a face of the VC.
# This happens if half of the lattice vector is outside the convex
# polytope defined by the rest of the vertices in `keep`.
bbt = basis @ basis.T
products = vertices @ bbt @ vertices.T
keep = np.array([True] * len(vertices))
for i in range(len(vertices)):
if not np.any(vertices[i]):
vertices[i] += 2 * np.array(displacements[i])
# Relevant if the projection of half of the lattice vector `vertices[i]`
# onto every other lattice vector in `veritces[keep]` is less than `0.5`.
mask = np.array([False if k == i else b for k, b in enumerate(keep)])
projections = 0.5 * products[i, mask] / np.diag(products)[mask]
if not np.all(np.abs(projections) < 0.5 - rtol):
keep[i] = False
vertices = vertices[keep]
vertices = np.concatenate([vertices, -vertices])
return vertices
......@@ -32,3 +32,35 @@ def test_cvp():
point = 50 * rng.randn(j)
assert np.array_equal(lll.cvp(point, mat, 10)[:3],
lll.cvp(point, mat, 3))
# Test equidistant vectors
# Cubic lattice
basis = np.eye(3)
vec = np.zeros((3))
assert len(lll.cvp(vec, basis, n=2, group_by_length=True)) == 7
assert len(lll.cvp(vec, basis, n=3, group_by_length=True)) == 19
vec = 0.5 * np.array([1, 1, 1])
assert len(lll.cvp(vec, basis, group_by_length=True)) == 8
vec = 0.5 * np.array([1, 1, 0])
assert len(lll.cvp(vec, basis, group_by_length=True)) == 4
vec = 0.5 * np.array([1, 0, 0])
assert len(lll.cvp(vec, basis, group_by_length=True)) == 2
# Square lattice with offset
offset = np.array([0, 0, 1])
basis = np.eye(3)[:2]
vec = np.zeros((3)) + rng.rand() * offset
assert len(lll.cvp(vec, basis, n=2, group_by_length=True)) == 5
assert len(lll.cvp(vec, basis, n=3, group_by_length=True)) == 9
vec = 0.5 * np.array([1, 1, 0]) + rng.rand() * offset
assert len(lll.cvp(vec, basis, group_by_length=True)) == 4
vec = 0.5 * np.array([1, 0, 0]) + rng.rand() * offset
assert len(lll.cvp(vec, basis, group_by_length=True)) == 2
# Hexagonal lattice
basis = np.array([[1, 0], [-0.5, 0.5 * np.sqrt(3)]])
vec = np.zeros((2))
assert len(lll.cvp(vec, basis, n=2, group_by_length=True)) == 7
assert len(lll.cvp(vec, basis, n=3, group_by_length=True)) == 13
vec = np.array([0.5, 0.5 / np.sqrt(3)])
assert len(lll.cvp(vec, basis, group_by_length=True)) == 3
assert len(lll.cvp(vec, basis, n=2, group_by_length=True)) == 6
assert len(lll.cvp(vec, basis, n=3, group_by_length=True)) == 12
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment