blaze2d 0.2.0

BLAZE - Band-structure LOBPCG Accelerated Zone Eigensolver for 2D Photonic Crystals

BLAZE - Band-structure LOBPCG Accelerated Zone Eigensolver

High-performance 2D photonic crystal band structure solver with Python bindings. Also supports Envelope Approximation (EA) eigenproblems for moiré lattice research.

Installation

pip install blaze2d

Quick Start

Maxwell Mode (Photonic Crystals)

from blaze import BulkDriver

# Load configuration from TOML file
driver = BulkDriver("config.toml")
print(f"Will run {driver.job_count} jobs, solver: {driver.solver_type}")

# Stream results as they complete (for live plotting)
for result in driver.run_streaming():
    if result['result_type'] == 'maxwell':
        print(f"Job {result['job_index']}: {result['num_bands']} bands")
        # result['bands'] is a 2D array [k_index][band_index]
        # result['distances'] is the k-path distance for plotting

# Or collect all results at once
results, stats = driver.run_collect()
print(f"Completed {len(results)} jobs in {stats['total_time_secs']:.2f}s")

EA Mode (Envelope Approximation for Moiré Lattices)

from blaze import BulkDriver

# Load EA configuration
driver = BulkDriver("ea_moire.toml")
print(f"Solver type: {driver.solver_type}")  # "ea"

# Stream eigenvalue results
for result in driver.run_streaming():
    if result['result_type'] == 'ea':
        eigenvalues = result['eigenvalues']
        print(f"Job {result['job_index']}: {len(eigenvalues)} eigenvalues")
        print(f"  Lowest 3: {eigenvalues[:3]}")
        print(f"  Converged: {result['converged']}")

Configuration

BLAZE uses TOML configuration files.

Maxwell Mode Example (photonic crystal sweep)

[bulk]
output_dir = "results"
output_mode = "full"

[bulk.ranges]
radius = { min = 0.2, max = 0.4, step = 0.05 }
eps_inside = { values = [9.0, 12.0, 13.0] }
polarization = ["te", "tm"]

[simulation]
resolution = 32
num_bands = 8

[lattice]
type = "square"

[[atoms]]
position = [0.5, 0.5]
radius = 0.3
epsilon = 1.0

EA Mode Example (moiré lattice)

[bulk]
threads = 4

[solver]
type = "ea"

[ea]
potential = "data/potential.bin"      # V(R) - binary f64 file
mass_inv = "data/mass_inv.bin"        # M⁻¹(R) - binary f64 file
# vg = "data/group_velocity.bin"      # Optional drift term
eta = 0.1
domain_size = [10.0, 10.0]

[grid]
nx = 64
ny = 64
lx = 10.0
ly = 10.0

[eigensolver]
n_bands = 12
tol = 1e-6
max_iter = 500

API Reference

BulkDriver

driver = BulkDriver(config_path: str, threads: int = 0)

• config_path: Path to TOML configuration file
• threads: Thread count (-1 = adaptive, 0 = auto, >0 = fixed)

Properties

Property	Type	Description
job_count	int	Number of jobs to execute
config_path	str	Path to configuration file
solver_type	str	"maxwell" or "ea"
is_ea	bool	True if EA solver

Methods

Method	Description
run_streaming()	Iterator yielding results as they complete
run_streaming_filtered(k_indices, band_indices)	Filtered streaming (Maxwell only)
run_collect(k_indices, band_indices)	Collect all results into list
run_batched(buffer_size_mb)	Run with buffered disk I/O
run()	Synchronous execution
dry_run()	Count jobs without executing

Result Dictionary

Each result is a dictionary. Common fields:

Key	Type	Description
job_index	int	Job identifier
result_type	str	"maxwell" or "ea"
params	dict	Simulation parameters
num_bands	int	Number of bands/eigenvalues

Maxwell-specific fields (result_type == "maxwell")

Key	Type	Description
k_path	list	K-points as (kx, ky) tuples
distances	list	Cumulative k-path distance
bands	list	2D array [k_index][band_index]
num_k_points	int	Number of k-points

EA-specific fields (result_type == "ea")

Key	Type	Description
eigenvalues	list	Computed eigenvalues (sorted)
n_iterations	int	LOBPCG iterations used
converged	bool	Whether solver converged
num_eigenvalues	int	Number of eigenvalues

Plotting Example

Maxwell Band Structure

import matplotlib.pyplot as plt
from blaze import BulkDriver

driver = BulkDriver("sweep.toml")

# Live plotting with streaming
plt.ion()
fig, ax = plt.subplots()

for result in driver.run_streaming():
    if result['result_type'] != 'maxwell':
        continue
    ax.clear()
    for band_idx in range(result['num_bands']):
        band = [result['bands'][k][band_idx] for k in range(result['num_k_points'])]
        ax.plot(result['distances'], band)
    ax.set_xlabel('k-path')
    ax.set_ylabel('ω (normalized)')
    ax.set_title(f"Job {result['job_index']}")
    plt.pause(0.01)

plt.ioff()
plt.show()

EA Eigenvalue Spectrum

import matplotlib.pyplot as plt
from blaze import BulkDriver

driver = BulkDriver("ea_config.toml")

results, stats = driver.run_collect()

# Plot eigenvalue spectrum
for result in results:
    eigenvalues = result['eigenvalues']
    plt.plot(eigenvalues, 'o-', label=f"Job {result['job_index']}")

plt.xlabel('Eigenvalue index')
plt.ylabel('Energy')
plt.title('EA Eigenvalue Spectrum')
plt.legend()
plt.show()

Filtered Streaming

For large simulations, filter to only the k-points and bands you need:

# Stream only Gamma (0), X (10), M (15) and first 4 bands
# Note: Filtering only applies to Maxwell results
for result in driver.run_streaming_filtered(
    k_indices=[0, 10, 15],
    band_indices=[0, 1, 2, 3]
):
    if result['result_type'] == 'maxwell':
        assert result['num_k_points'] == 3
        assert result['num_bands'] == 4

EA Input Data Format

EA mode requires binary input files for the effective Hamiltonian components. All files are raw f64 (little-endian) in row-major (C-order):

• Potential (potential.bin): Nx × Ny values = V(x, y)
• Mass inverse (mass_inv.bin): Nx × Ny × 4 values = [m_xx, m_xy, m_yx, m_yy]
• Group velocity (vg.bin, optional): Nx × Ny × 2 values = [vg_x, vg_y]

Python example to generate input data:

import numpy as np

nx, ny = 64, 64
Lx, Ly = 10.0, 10.0

# Coordinate grids
x = np.linspace(0, Lx, nx, endpoint=False)
y = np.linspace(0, Ly, ny, endpoint=False)
X, Y = np.meshgrid(x, y, indexing='ij')

# Periodic potential
G = 2 * np.pi / Lx
V = 0.1 * (np.cos(G * X) + np.cos(G * Y))

# Isotropic mass (m* = 0.05)
m_star = 0.05
m_inv = np.zeros((nx, ny, 4))
m_inv[:, :, 0] = 1.0 / m_star  # m_xx
m_inv[:, :, 3] = 1.0 / m_star  # m_yy

# Save
V.astype(np.float64).tofile('data/potential.bin')
m_inv.astype(np.float64).tofile('data/mass_inv.bin')

License

MIT