Simple Interpolation
One of the main purposes of brille is to perform linear interpolation. This short tutorial will teach you how to perform simple linear interpolation using brille. All examples here are used in the continuous testing of brille and are from the file wrap/tests/test_4_interpolation.py available in the repository.
Interpolation of a scalar field
A scalar field is any function \(\phi(\mathbf{Q})\) which is single-valued at all \(\mathbf{Q}\). Here we will interpolate a scalar field of the form \(\phi(\mathbf{Q}) = \mathbf{Q}\cdot\hat{\mathbf{x}}\), which is a linear scalar field.
Construct the lattice
Linear interpolation can only be performed using brille if a grid is constructed first. The construction of a grid requires that we define a space-filling lattice, so we start with that:
from brille import Lattice
lattice = Lattice((1, 1, 1), (90, 90, 90), real_space=False)
The created lattice spans reciprocal space with three orthogonal basis vectors each of unit length, \(1 \AA^{-1}\), and as a result has a unit cell with volume \(1 \AA^{-3}\) which we can access through lattice.volume.
Construct the Brillouin zone polyhedron
To construct the grid we must define its extent, which in brille is always a
first or irreducible Brillouin zone polyhedron. Finding this polyhedron may not
be trivial, so it is an independent step (with options to control the algorithm
behaviour, see BrillouinZone for details).
from brille import BrillouinZone
brillouin_zone = BrillouinZone(lattice)
We did not specify any symmetry information for this lattice, so the default \(P 1\) is used, which has an identical first and irreducible Brillouin zone. Therefore brillouin_zone contains an ‘irreducible’ polyhedron which identical to the first Brillouin zone polyhedron, accessible through brillouin_zone.ir_polyhedron and brillouin_zone.polyhedron, respectively.
Construct the interpolation grid
Finally, we can construct the interpolation grid. We have multiple options
within brille, but nearly always we should use one of the Hybrid grids.
In this short tutorial we will interpolate purely real-valued data, so we pick
a BZTrellisQdd which holds real-valued eigenvalues
and eigenvectors for every vertex in its grid.
To create a grid we must specify an upper-limit to the grid-cell size, which
also specifies a lower-limit to the number of points in the grid.
Each grid cell contributes one or more vertex to the grid, with the exact
number dependent on details of the grid and bounding polyhedron.
Here we will aim for at least \(100\) grid points, or a maximum cell
volume of \(1\times10^{-2} \AA^{-3}\).
from brille import BZTrellisQdd
maximum_cell_volume = lattice.volume / 100
grid = BZTrellisQdd(brillouin_zone, maximum_cell_volume)
The points in the grid are accessible through grid.rlu or grid.invA in units of the reciprocal lattice units or \(\AA^{-1}\), respectively; in this case the two units are the same.
Provide the data to be interpolated
The grid must be provided with data corresponding to each grid point in order
to be used to perform linear interpolation. This is accomplished using one of
the fill() method overloads.
Note
The fill methods require input for both the eigenvalues and eigenvectors
to be interpolated. This specialization is a limitation which may be removed
in a future version of brille.
Here we will provide the same data as eigenvalue and eigenvector. In addition we must provide metadata for the eigenvalues and eigenvectors stored at each grid point. The metadata is a list-like object consisting of the number of scalar values, number of vector-elements, and number of matrix-elements stored at each grid point for each eigenvalue-eigenvector pair (if a system has multiple branches there will be multiple pairs). A fourth piece of metadata determines how the vector- and matrix- elements are handled when symmetry operations are used in the interpolation. The scalar function we want to interpolate has one scalar value per grid point, and we will not use symmetry operations in the interpolation so the default, \(0\), is acceptable here.
def phi(q):
return q[:, 0]
phi_of_q = phi(grid.invA)
metadata = (1, 0, 0, 0)
grid.fill(phi_of_q, metadata, phi_of_q, metadata)
Perform the interpolation
The extent of the grid is the first Brillouin zone polyhedron which, in this case, has all coordinates \(x_i \in (-0.5, 0.5)\). If the grid is used to interpolate at a \(\mathbf{Q}\) inside of the first Brillouin zone, the result should be identical to \(\phi(\mathbf{Q})\); but if \(\mathbf{Q}\) is not in the first Brillouin zone the result will be \(\phi(\mathbf{Q} - \mathbf{G})\) where \(\mathbf{G}\) is the nearest integer lattice point to \(\mathbf{Q}\).
import numpy as np
q_pts = np.random.rand(10,3) - 0.5 # all in range (-0.5, 0.5)
interpolated_values, interpolated_vectors = grid.interpolate_at(q_pts)
if np.allclose(interpolated_values, interpolated_vectors):
print('The interpolated values and vectors match!')
else:
print('This should be an error, and impossible.')
if np.allclose(interpolated_values, phi(q_pts)):
print('The interpolation worked as expected!')
else:
print('Linear interpolation of a linear scalar field did not work?!')
q_pts *= 20 # now all in the range (-10, 10)
interpolated_values, interpolated_vectors = grid.interpolate_at(q_pts)
if np.allclose(interpolated_values, phi(q_pts - np.round(q_pts))):
print('The interpolation still worked as expected!')
else:
print('Linear interpolation or subtracting G did not work?!')
The interpolated results should be exact (to machine precision) in this case since the linear interpolation of a linear function is exactly independent of the interpolation grid step size. Interpolating any non-linear function will naturally introduce some error in its estimate of the function.
Interpolation of a dispersing excitation
Excitations which have energies dependent on \(\mathbf{Q}\) are said to be dispersive. One such dispersive excitation is the accoustic ferromagnetic spinwave in iron, which has a dispersion
where \(\delta\) is a single-ion anisotropy term, \(J\) is the exchange energy, \(a\) is the iron lattice parameter, and the \(\hat{\mathbf{x}}_i\) are the three Cartesian directions.
The spinwave dispersion has the same periodicity as the iron crystal lattice so
it is possible to use brille to estimate its energy at arbitrary
\(\mathbf{Q}\) using linear interpolation.
Setup
The lattice parameter of iron is \(a = 2.87 \AA\) and its cubic lattice has the symmetry of the spacegroup with Hermann-Mauguin symbol \(I m \bar{3} m\). We will construct an irreducible Brillouin zone polyhedron, and then use it to produce a hybrid interpolation grid with at least \(1000\) points.
from brille import Lattice, BrillouinZone, BZTrellisQdd
a_fe = 2.87
direct_lattice = Lattice((a_fe, a_fe, a_fe), (90, 90, 90), 'I m -3 m')
brillouin_zone = BrillouinZone(direct_lattice.star)
grid = BZTrellisQdd(brillouin_zone, brillouin_zone.ir_polyhedron.volume/1000)
Fill
We then can use the relative lattice unit grid points to fill the dispersion information into the grid. Again we will provide the same information for the eigenvalues and eigenvectors, and again there is a single scalar value for the single dispersion branch. The default fourth metadata value would not be appropriate for interpolation of true eigenvectors in this case, but has no effect for scalars and so does not impact this case.
def omega_fe(q, J=16, delta=0.01):
from numpy import cos, pi, prod
return delta + 8*J*(1 - prod(cos(pi*q), axis=1))
omega = omega_fe(grid.rlu)
metadata = (1, 0, 0, 0)
grid.fill(omega, metadata, omega, metadata)
Interpolate
Since we have provided symmetry information and produced a grid bounded by
an irreducible Brillouin zone polyhedron, we must use the irreducible
interpolation method of the grid, ir_interpolate_at().
As the dispersion relation is not a linear function, the linear interpolation will always introduce error when compared to the true function. The size of this error depends on the grid vertex spacing and the local curvature of the dispersion. Here we can only verify that the interpolation returns the stored values at the grid points.
from numpy import allclose
values, vectors = grid.ir_interpolate_at(grid.rlu)
if not allclose(values, omega) or not allclose(values, grid.values):
print("This should be an error")
if not allclose(vectors, omega) or not allclose(vectors, grid.vectors):
print("This should be an error")
With the assurance that the interpolation works we can examine more interesting paths through reciprocal space:
from matplotlib.pyplot import plot, gca, setp
from numpy import linspace, vstack, array, arange
x = array([[0, 0, 0], [1, 0, 0], [0.5, 0.5, 0], [0, 0, 0], [0.5, 0.5, 0.5]])
n_pts = 100
path = vstack([linspace(x[i], x[i+1], n_pts) for i in range(len(x)-1)])
values, _ = grid.ir_interpolate_at(path)
x_plot = arange(len(path))
ticks_at = n_pts * arange(len(x))
tick_labels = [f'({z[0]} {z[1]} {z[2]})' for z in x]
plot(x_plot, values, '-k', x_plot, omega_fe(path), '--r')
setp(gca(), xticks=ticks_at, xticklabels=tick_labels)
setp(gca(), ylabel=r'$\omega(\mathbf{Q})$', xlabel=r'$\mathbf{Q}$')
Produces the following plot:
Interpolated iron spinwave dispersion (black) and exact solution (red).
On closer inspection, the linear-interpolation introduced errors are more obvious:
plot(x_plot, values, '-k', x_plot, omega_fe(path), '--r')
setp(gca(), xticks=ticks_at, xticklabels=tick_labels)
setp(gca(), ylabel=r'$\omega(\mathbf{Q})$', xlabel=r'$\mathbf{Q}$')
setp(gca(), xlim=(80, 120), ylim=(247,257))
Enlarged region of interpolated iron spinwave dispersion (black) and exact solution (red).