Switch from BZTrellisQ to BZMeshQ¶

BZMeshQdc and its siblings can replace BZTrellisQdc in most code: they are filled and interpolated in the same way. This page shows what to change, what to watch for, and how to check whether the mesh suits your system.

The code on this page comes from trellis_to_mesh.py, which the tests run.

Change the constructor¶

A trellis like this:

volume = bz.ir_polyhedron.volume
trellis = BZTrellisQdc(bz, node_volume_fraction=volume / points)
trellis.fill(*grid_data(trellis.rlu))

becomes:

mesh = BZMeshQdc(bz, max_size=volume / (6 * points))
mesh.fill(*grid_data(mesh.rlu))

Both size parameters are volumes in Å⁻³, not fractions:

• The trellis's node_volume_fraction is, despite its name, the volume of one cube of the trellis.
• The mesh's max_size is the largest volume of one tetrahedron.

A cell of the mesh's grid holds six tetrahedra, so dividing by six gives a mesh with about as many vertices as the trellis:

Vertices for node_volume_fraction = V/N and max_size = V/(6N), where V is the irreducible zone's volume

All other arguments, and fill, ir_interpolate_at, sort, and saving and loading, work as before.

If you create grids with brille.utils.create_grid, the mesh is its default from brille 0.9: pass max_size in place of node_volume_fraction. Passing node_volume_fraction, or trellis=True, still makes a trellis.

grid = create_grid(bz, complex_vectors=True, max_size=volume / (6 * points))

Through Euphonic¶

Euphonic's BrilleInterpolator.from_force_constants takes grid_type="mesh". For the same grid_npts, it asks the mesh for between a fifth and a half as many points as the trellis (in Euphonic 2.1). To keep the number of points, give the mesh's size directly:

from brille import BrillouinZone
from euphonic.brille import BrilleInterpolator

volume = BrillouinZone(lattice).ir_polyhedron.volume   # your crystal's lattice
interpolator = BrilleInterpolator.from_force_constants(
    force_constants, grid_type="mesh", grid_kwargs={"max_size": volume / (6 * grid_npts)})

What else changes¶

The mesh adds:

• refinement_points and refine, to add vertices where the interpolation needs them (see Building and refining a mesh);
• max_points, which caps the number of vertices.

Interpolated values change slightly, because the vertices are different: by about the interpolation error, not more. On the zone boundary, the mesh matches itself across equivalent faces, so interpolation is continuous there.

A saved trellis can only be read back as a trellis, so rebuild your grids as meshes rather than converting files. Meshes saved by earlier versions of brille, which built them with TetGen, still load, but cannot be refined.

Trade-offs¶

These figures compare the two grids for three lattices. They were measured on one machine, with one thread, interpolating a smooth symmetric function at 100,000 random points. Accuracy is given as rms error × vertices^2/3: for linear interpolation of a smooth function this hardly depends on the number of vertices, so it compares grids of different sizes, and lower is better.

In short:

• Interpolation with the mesh is 2–3 times faster: it finds the tetrahedron holding a point in constant time.
• Building the mesh is faster, by a factor that grows with its size: 17 times for the triclinic lattice above.
• Accuracy per vertex is the same for fine grids, or better for the mesh (fcc above). For coarse grids of low-symmetry lattices, a trellis of a few hundred vertices can be up to twice as accurate per vertex as a mesh of the same size.
• Memory held by the mesh is similar to the trellis's, or less for the same number of vertices.

Stay with the trellis if you rely on its node methods or on interpolate_at, if you need results identical to earlier runs, or if you use very coarse grids of low-symmetry lattices. Otherwise, the mesh is faster, and refinement lets it put vertices where your model needs them.

Measure for your own system¶

The best test is your own model. This builds both grids at three sizes, then reports build time, interpolation time and rms error; replace lattice and model in the script with yours.

def measure(make, q, repeat=3):
    """Build a grid, fill it with the model, and time interpolating at q"""
    start = time.perf_counter()
    grid = make()
    built = time.perf_counter() - start
    grid.fill(*grid_data(grid.rlu))
    best = np.inf
    for _ in range(repeat):
        start = time.perf_counter()
        values = grid.ir_interpolate_at(q, threads=1)[0][:, 0]
        best = min(best, time.perf_counter() - start)
    error = np.sqrt(np.mean((values - model(q)) ** 2))
    return len(grid.rlu), built, 1e6 * best / len(q), error


def compare(sizes=(1000, 5000, 20000)):
    """Build a trellis and a mesh of about each size in sizes, and compare them"""
    q = np.random.default_rng(1).uniform(-1, 1, (20000, 3))
    print(f"{'grid':8s} {'vertices':>8s} {'build/s':>8s} {'µs/point':>9s} {'rms error':>10s}")
    for points in sizes:
        for name, make in (("trellis", lambda: BZTrellisQdc(bz, node_volume_fraction=volume / points)),
                           ("mesh", lambda: BZMeshQdc(bz, max_size=volume / (6 * points)))):
            vertices, built, per_point, error = measure(make, q)
            print(f"{name:8s} {vertices:8d} {built:8.2f} {per_point:9.2f} {error:10.2e}")


if __name__ == "__main__":
    compare()

For the example's triclinic lattice:

grid     vertices  build/s  µs/point  rms error
trellis      1941     0.71      6.19   2.61e-02
mesh         2886     0.48      1.55   2.53e-02
trellis      7787     3.10      6.61   9.11e-03
mesh         9947     1.12      2.15   9.06e-03
trellis     26432    12.76      7.30   3.86e-03
mesh        32505     2.76      3.11   3.59e-03