blaze2d 0.3.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.

v0.3.0: EA mode now emits eigenvectors alongside eigenvalues!

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
import numpy as np

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

# Stream eigenvalue and eigenvector results
for result in driver.run_streaming():
    if result['result_type'] == 'ea':
        eigenvalues = result['eigenvalues']
        eigenvectors = result['eigenvectors']  # List of complex 2D arrays
        grid_dims = result['grid_dims']  # [nx, ny]
        
        print(f"Job {result['job_index']}: {len(eigenvalues)} eigenvalues")
        print(f"  Lowest 3: {eigenvalues[:3]}")
        print(f"  Converged: {result['converged']}")
        
        # Eigenvectors are complex arrays of shape (nx*ny, 2) -> [re, im]
        # Reshape to 2D field:
        psi_0 = np.array(eigenvectors[0])  # First eigenvector
        psi_0_complex = psi_0[:, 0] + 1j * psi_0[:, 1]
        psi_0_2d = psi_0_complex.reshape(grid_dims)

Configuration

BLAZE uses TOML configuration files.

Maxwell Mode Example (photonic crystal sweep)

# Section-less entries at top
polarization = "TM"

[bulk]
threads = 4
verbose = true

[grid]
nx = 32
ny = 32
lx = 1.0
ly = 1.0

[geometry]
eps_bg = 12.0

[geometry.lattice]
type = "square"
a = 1.0

[[geometry.atoms]]
pos = [0.5, 0.5]
radius = 0.3
eps_inside = 1.0

[path]
preset = "square"
segments_per_leg = 12

[eigensolver]
n_bands = 8
tol = 1e-6
max_iter = 200

[ranges]
eps_bg = { min = 10.0, max = 14.0, step = 1.0 }
polarization = ["TM", "TE"]

[output]
mode = "full"
directory = "./results"
prefix = "job"

EA Mode Example (moiré lattice)

# EA mode requires [solver] type = "ea"
# Eigenvectors are always emitted alongside eigenvalues

[bulk]
threads = 4
verbose = true

[solver]
type = "ea"

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

[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]
periodic = true

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

[output]
mode = "full"
directory = "./ea_output"
prefix = "ea_job"

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)
eigenvectors	list	List of eigenvectors, each as (N, 2) array [re, im]
grid_dims	list	Grid dimensions [nx, ny] for reshaping eigenvectors
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 and Eigenvector Visualization

import matplotlib.pyplot as plt
import numpy as np
from blaze import BulkDriver

driver = BulkDriver("ea_config.toml")
results, stats = driver.run_collect()

result = results[0]
eigenvalues = result['eigenvalues']
eigenvectors = result['eigenvectors']
nx, ny = result['grid_dims']

# Plot eigenvalue spectrum
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].plot(eigenvalues, 'o-')
axes[0].set_xlabel('Eigenvalue index')
axes[0].set_ylabel('Energy')
axes[0].set_title('EA Eigenvalue Spectrum')

# Plot ground state |ψ|²
psi_0 = np.array(eigenvectors[0])  # Shape: (nx*ny, 2)
psi_0_complex = psi_0[:, 0] + 1j * psi_0[:, 1]
psi_0_2d = psi_0_complex.reshape((nx, ny))

im = axes[1].imshow(np.abs(psi_0_2d)**2, origin='lower', cmap='hot')
axes[1].set_title(f'Ground state |ψ₀|² (E = {eigenvalues[0]:.4f})')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
plt.colorbar(im, ax=axes[1])

plt.tight_layout()
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