Usage

Here we demonstrate the basic usage for the SmoQyKPMCore package. Let us first import the relevant packages we will want to use in this example.

using LinearAlgebra
using SparseArrays

Package for making figures

using CairoMakie
CairoMakie.activate!(type = "svg")

using SmoQyKPMCore

In this usage example we will consider a 1D chain tight-binding model

\[\hat{H} = -t \sum_{i,\sigma} (\hat{c}_{i+1,\sigma}^\dagger \hat{c}_{i,\sigma} + {\rm H.c.}) = \sum_\sigma \hat{c}_{i,\sigma}^\dagger [H_{i,j}] \hat{c}_{j,\sigma},\]

where $\hat{c}_{i,\sigma}^\dagger \ (\hat{c}_{i,\sigma})$ creates (annihilates) a spin-$\sigma$ electron on site $i$ in the lattice, and $t$ is the nearest-neighbor hopping amplitude.

Density Matrix Approximation

With $H$ the Hamiltonian matrix, the corresponding density matrix is given by $\rho = f(H)$ where

\[f(\epsilon) = \frac{1}{1+e^{\beta (\epsilon-\mu)}} = \frac{1}{2} \left[ 1 + \tanh\left(\tfrac{\beta(\epsilon-\mu)}{2}\right) \right]\]

is the Fermi function, with $\beta = 1/T$ the inverse temperature and $\mu$ the chemical potential. Here we will applying the kernel polynomial method (KPM) algorithm to approximate the density matrix $\rho$ by a Chebyshev polynomial expansion.

Let us first define the relevant model parameter values.

# Nearest-neighbor hopping amplitude.
t = 1.0

# Chemical potential.
μ = 0.0

# System size.
L = 16

# Inverse temperature.
β = 4.0;

Now we need to construct the Hamiltonian matrix $H$.

H = zeros(L,L)
for i in 1:L
    j = mod1(i+1,L)
    H[j,i], H[i,j] = -t, -t
end
H

16×16 Matrix{Float64}:
  0.0  -1.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0  -1.0
 -1.0   0.0  -1.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0  -1.0   0.0  -1.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0  -1.0   0.0  -1.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0  -1.0   0.0  -1.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0  -1.0   0.0  …   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0  -1.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0     -1.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …   0.0  -1.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0     -1.0   0.0  -1.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0  -1.0   0.0  -1.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0  -1.0   0.0  -1.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0  -1.0   0.0  -1.0
 -1.0   0.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0  -1.0   0.0

We also need to define the Fermi function.

# Fermi function.
fermi(ϵ, μ, β) = (1+tanh(β*(ϵ-μ)/2))/2

fermi (generic function with 1 method)

Let us calculate the exact density matrix $\rho$ so that we may asses the accuracy of the KPM method.

# Diagonalize the Hamiltonian matrix.
ϵ, U = eigen(H)

# Calculate the density matrix.
ρ = U * Diagonal(fermi.(ϵ, μ, β)) * adjoint(U)

16×16 Matrix{Float64}:
  0.5          -0.30913      -2.70673e-16  …   4.16315e-16  -0.30913
 -0.30913       0.5          -0.30913          0.0826932     2.56734e-16
 -2.42894e-16  -0.30913       0.5              6.92028e-18   0.0826932
  0.0826932    -2.91452e-16  -0.30913         -0.0320635     2.77509e-17
  1.17929e-16   0.0826932    -3.26184e-16      1.38592e-17  -0.0320635
 -0.0320635     9.0187e-17    0.0826932    …   0.00866812   -9.71492e-17
 -1.73505e-16  -0.0320635     2.21989e-16      1.80393e-16   0.00866812
  0.00866812   -1.04102e-16  -0.0320635        0.00866812    6.93843e-17
  1.38745e-16   0.00866812   -1.31895e-16     -1.52674e-16   0.00866812
  0.00866812    2.01209e-16   0.00866812      -0.0320635    -1.17966e-16
 -1.39104e-17   0.00866812   -2.08223e-16  …   1.3182e-16   -0.0320635
 -0.0320635     6.24314e-17   0.00866812       0.0826932     2.0812e-17
  1.38452e-17  -0.0320635     2.08111e-16     -2.63697e-16   0.0826932
  0.0826932    -7.21664e-16  -0.0320635       -0.30913       2.8449e-16
  4.16301e-16   0.0826932     2.07608e-17      0.5          -0.30913
 -0.30913       2.49782e-16   0.0826932    …  -0.30913       0.5

Now let use the KPMExpansion type to approximate $\rho$. We will need to define the order $M$ of the expansion and give approximate bounds for the eigenspectrum of $H$, making sure to overestimate the the true interval spanned by the eigenvalues of $H$. Because we are considering the simple non-interacting model, the exact eigenspectrum of $H$ is known and spans the interval $\epsilon \in [-2t, 2t].$

# Define eigenspectrum bounds.
bounds = (-3.0t, 3.0t)

# Define order of Chebyshev expansion used in KPM approximation.
M = 10

# Initialize expansion.
kpm_expansion = KPMExpansion(x -> fermi(x, μ, β), bounds, M);

Now let us test and see how good a job it does approximating the density matrix $\rho$, using the kpm_eval! function to do so.

# Initialize KPM density matrix approxiatmion
ρ_kpm = similar(H)

# Initialize additional matrices to avoid dynamic memory allocation.
mtmp = zeros(L,L,3)

# Calculate KPM density matrix approximation.
kpm_eval!(ρ_kpm, H, kpm_expansion, mtmp)

# Check how good an approximation it is.
println("Matrix Error = ", norm(ρ_kpm - ρ)/norm(ρ) )

Matrix Error = 0.01669022782185739

It is also possible to efficiently multiply vectors by KPM approximation to the density matrix $\rho$. Let us test this functionality with the kpm_mul! on a random vector, seeing how accurate the result is.

# Initialize random vector.
v = randn(L)

# Calculate exact product with density matrix.
ρv = ρ * v

# Initialize vector to contain approximate product with densith matrix.
ρv_kpm = similar(v)

# Initialize to avoid dynamic memory allocation.
vtmp = zeros(L, 3)

# Calculate approximate product with density matrix.
kpm_mul!(ρv_kpm, H, kpm_expansion, v, vtmp)

# Check how good the approximation is.
println("Vector Error = ", norm(ρv_kpm - ρv) / norm(ρv) )

Vector Error = 0.020833304968523727

Let us now provide the KPM approximation with better bounds on the eigenspectrum and also increase the order the of the expansion and see how the result improves. We will use the kpm_update! function to update the KPMExpansion in-place.

# Define eigenspectrum bounds.
bounds = (-2.5t, 2.5t)

# Define order of Chebyshev expansion used in KPM approximation.
M = 100

# Update KPM approximation.
kpm_update!(kpm_expansion, x -> fermi(x, μ, β), bounds, M)

# Calculate KPM density matrix approximation.
kpm_eval!(ρ_kpm, H, kpm_expansion, mtmp)

# Check how good an approximation it is.
println("Matrix Error = ", norm(ρ_kpm - ρ)/norm(ρ) )

# Calculate approximate product with density matrix.
kpm_mul!(ρv_kpm, H, kpm_expansion, v, vtmp)

# Check how good the approximation is.
println("Vector Error = ", norm(ρv_kpm - ρv) / norm(ρv) )

Matrix Error = 1.031960839804449e-14
Vector Error = 7.336124683100945e-15

Now let me quickly demonstrate how we may approximate the the trace $\rho$ using a set of random vector $R_n$ using the kpm_dot function.

# Calculate exact trace.
trρ = tr(ρ)
println("Exact trace = ", trρ)

# Number of random vectors
N = 100

# Initialize random vectors.
R = randn(L, N)

# Initialize array to avoid dynamic memory allocation
Rtmp = zeros(L, N, 3)

# Approximate trace of density matrix.
trρ_approx = kpm_dot(H, kpm_expansion, R, Rtmp)
println("Approximate trace = ", trρ_approx)

# Report the error in the approximation.
println("Trace estimate error = ", abs(trρ_approx - trρ)/trρ)

Exact trace = 8.000000000000002
Approximate trace = 8.19195416031251
Trace estimate error = 0.02399427003906362

Density of States Approximation

Next let us demonstrate how we may approximate the density of states $\mathcal{N}(\epsilon)$ for a 1D chain tight-binding model. The first step a very larger Hamiltonian matrix $H,$ which we represent as a spare matrix.

# Size of 1D chain considered
L = 10_000

# Construct sparse Hamiltonian matrix.
rows = Int[]
cols = Int[]
for i in 1:L
    j = mod1(i+1, L)
    push!(rows,i)
    push!(cols,j)
    push!(rows,j)
    push!(cols,i)
end
vals = fill(-t, length(rows))
H = sparse(rows, cols, vals)

10000×10000 SparseArrays.SparseMatrixCSC{Float64, Int64} with 20000 stored entries:
⎡⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⎤
⎢⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥
⎣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦

Let us now calculate the fist $M$ moments $\mu_m$ using $N$ random vectors using the function kpm_moments.

# Number of moments to calculate.
M = 250

# Number of random vectors used to approximate moments.
N = 100

# Initialize random vectors.
R = randn(L, N)

# Calculate the moments.
μ_kpm = kpm_moments(M, H, bounds, R);

Having calculate the moments, let us next evaluate the density of states $\mathcal{N}(\epsilon)$ at $P$ points with the kpm_density function.

# Number of points at which to evaluate the density of states.
P = 1000

# Evaluate density of states.
dos, ϵ = kpm_density(P, μ_kpm, bounds);

Without regularization, the approximate for the density of states generated above will have clearly visible Gibbs oscillations. To partially suppress these artifacts, we apply the Jackson kernel to the moments $\mu$ using the apply_jackson_kernel function.

# Apply Jackson kernel.
μ_jackson = apply_jackson_kernel(μ_kpm)

# Evaluate density of states.
dos_jackson, ϵ_jackson = kpm_density(P, μ_jackson, bounds);

Having approximated the density of states, let us now plot it.

fig = Figure(
    size = (700, 400),
    fonts = (; regular= "CMU Serif"),
    figure_padding = 10
)

ax = Axis(
    fig[1, 1],
    aspect = 7/4,
    xlabel = L"\epsilon", ylabel = L"\mathcal{N}(\epsilon)",
    xticks = range(start = bounds[1], stop = bounds[2], length = 5),
    xlabelsize = 30, ylabelsize = 30,
    xticklabelsize = 24, yticklabelsize = 24,
)

lines!(
    ϵ, dos,
    linewidth = 2.0,
    alpha = 1.0,
    color = :red,
    linestyle = :solid,
    label = "Dirichlet"
)

lines!(
    ϵ_jackson, dos_jackson,
    linewidth = 3.0,
    alpha = 1.0,
    color = :black,
    linestyle = :solid,
    label = "Jackson"
)

xlims!(
    ax, bounds[1], bounds[2]
)

ylims!(
    ax, 0.0, 1.05 * maximum(dos)
)

axislegend(
    ax, halign = :center, valign = :top, labelsize = 30
)

fig