# Examples

The following are several examples of how to use the algorithms and integrands provided by AutoBZCore.jl. For background on the essential interface see the Problem definitions page

## Green's function integration

A common integral appearing in Dynamical mean-field theory is that of the local Green's function:

\[G(\omega) = \int d^d \vec{k}\ \operatorname{Tr} \left[ \left( \omega - H \left( \vec{k} \right) - \Sigma(\omega) \right)^{-1} \right].\]

For simplicity, we take $\Sigma(\omega) = -i\eta$. We can define the integrand as a function of $\vec{k}$ and $H$ and parameters $\eta, \omega$.

```
using LinearAlgebra
gloc_integrand(k, (; h, η, ω)) = tr(inv(complex(ω,η)*I-h(k)))
```

`gloc_integrand (generic function with 1 method)`

Here we use named tuple destructuring syntax to unpack a named tuple of parameters in the function definition.

Commonly, $H(\vec{k})$ is evaluated using Wannier interpolation, i.e. as a Fourier series. For a simple tight-binding model, the integer lattice, the Hamiltonian is given by

\[H(k) = \cos(2\pi k) = \frac{1}{2} \left( e^{2\pi ik} + e^{-2\pi ik} \right)\]

We can use the built-in function `cos`

to evaluate this, however, for more complex Fourier series it becomes easier to use the representation in terms of Fourier coefficients. Using the package FourierSeriesEvaluators.jl, we can define $H(k) = \cos(2\pi k)$ by the following:

```
using FourierSeriesEvaluators
h = FourierSeries([0.5, 0.0, 0.5]; period=1, offset=-2)
```

`3-element and (1,)-periodic FourierSeries with Float64 coefficients, (0,) derivative, (-2,) offset`

The coefficient values of $1/2$ can be determined from Euler's formula, as used in the expansion of $\cos$ above, and the value of `offset`

is chosen to offset the coefficient array indices, `1:3`

since Julia has 1-based indexing, to correspond to values of $n$ in the phase factors $e^{2\pi i n k}$ used in the Fourier series above, i.e. `-1:1`

. Now we proceed to the define the `IntegralProblem`

and solve it with a generic adaptive integration scheme, `QuadGKJL`

```
using AutoBZCore
dom = (0, 1)
p = (; h, η=0.1, ω=0.0)
prob = IntegralProblem(gloc_integrand, dom, p)
alg = QuadGKJL()
solve(prob, alg; abstol=1e-3).value
```

`4.327818095284375e-7 - 0.995037545189551im`

## BZ integration

To perform integration over a Brillouin zone, we can load one using the `load_bz`

function and then construct an `AutoBZProblem`

to solve. Since the Brillouin zone may be reduced using point group symmetries, a common optimization, it is also required to specify the symmetry representation of the integrand. Continuing the previous example, the trace of the Green's function has no band/orbital degrees of freedom and transforms trivially under the point group, so it is a `TrivialRep`

. The previous calculation can be replicated as:

```
using AutoBZCore
bz = load_bz(FBZ(), 2pi*I(1))
p = (; h, η=0.1, ω=0.0)
prob = AutoBZProblem(TrivialRep(), IntegralFunction(gloc_integrand), bz, p)
alg = TAI()
solve(prob, alg; abstol=1e-3).value
```

`5.551115123125783e-17 - 0.9950375451895511im`

Now we proceed to multi-dimensional integrals. In this case, Wannier interpolation is much more efficient when Fourier series are evaluated one variable at a time. To understand, this suppose we have a series defined by $M \times M$ coefficients (i.e. a 2d series) that we want to evaluate on an $N \times N$ grid. Naively evaluating the series at each grid point will require $\mathcal{O}(M^{2} N^{2})$ operations, however, we can reduce the complexity by pre-evaluating certain coefficients as follows

\[f(x, y) = \sum_{m,n=1}^{M} f_{nm} e^{i(nx + my)} = \sum_{n=1}^{M} e^{inx} \left( \sum_{m=1}^{M} f_{nm} e^{imy} \right) = \sum_{n=1}^{M} e^{inx} \tilde{f}_{n}(y)\]

This means we can evaluate the series on the grid in $\mathcal{O}(M N^2 + M^2 N)$ operations. When $N \gg M$, this is $\mathcal{O}(M N^{2})$ operations, which is comparable to the computational complexity of a multi-dimensional FFT. Since the constants of a FFT may not be trivial, this scheme is competitive.

Let's use this with a Fourier series corresponding to $H(\vec{k}) = \cos(2\pi k_{x}) + \cos(2\pi k_{y})$ and define a new method of `gloc_integrand`

that accepts the (efficiently) evaluated Fourier series in the second argument

```
h = FourierSeries([0.0; 0.5; 0.0;; 0.5; 0.0; 0.5;; 0.0; 0.5; 0.0]; period=1, offset=-2)
gloc_integrand(k, h_k, (; η, ω)) = tr(inv(complex(ω,η)*I-h_k))
```

`gloc_integrand (generic function with 2 methods)`

Similar to before, we construct an `AutoBZCore.IntegralProblem`

and this time we take the integration domain to correspond to the full Brillouin zone of a square lattice with lattice vectors `2pi*I(2)`

.

```
integrand = FourierIntegralFunction(gloc_integrand, h)
bz = load_bz(FBZ(2), 2pi*I(2))
p = (; η=0.1, ω=0.0)
prob = AutoBZProblem(TrivialRep(), integrand, bz, p)
alg = IAI()
solve(prob, alg).value
```

`1.3530843112619095e-16 - 1.3941704645339514im`

This package provides several `AutoBZProblem`

algorithms that we can use to solve the multidimensional integral.

The repo's demo on density of states provides a complete example of how to compute and interpolate an integral as a function of its parameters using the `init`

and `solve!`

interface

## Density of States

Computing the density of states (DOS) of a self-adjoint, or Hermitian, operator is a related, but distinct problem to the integrals also presented in this package. In fact, many DOS algorithms will compute integrals to approximate the DOS of an operator by introducing an artificial broadening. To handle the $T=0^{+}$ limit of the broadening, we implement the well-known Gilat-Raubenheimer method as an algorithm for the `AutoBZCore.DOSProblem`

Using the `AutoBZCore.init`

and `AutoBZCore.solve!`

functions, it is possible to construct a cache to solve a `DOSProblem`

for several energies or several Hamiltonians. As an example of solving for several energies,

```
using AutoBZCore, FourierSeriesEvaluators, StaticArrays
h = FourierSeries(SMatrix{1,1,Float64,1}.([0.5, 0.0, 0.5]), period=1.0, offset=-2)
E = 0.3
bz = load_bz(FBZ(), [2pi;;])
prob = DOSProblem(h, E, bz)
alg = GGR(; npt=100)
cache = init(prob, alg)
Es = range(-1, 1, length=10) * 1.1
data = map(Es) do E
cache.domain = E
solve!(cache).value
end
```

```
10-element Vector{Float64}:
0.0
0.5940536926606439
0.3934526572333864
0.34235097468581055
0.3208398038674902
0.32083980386749017
0.34235097468581066
0.3934526572333864
0.5940536926606439
0.0
```

As an example of interchanging Hamiltonians, which must remain the same type, we can double the energies, which will halve the DOS

```
cache.domain = E
sol1 = AutoBZCore.solve!(cache)
h.c .*= 2
cache.isfresh = true
cache.domain = 2E
sol2 = AutoBZCore.solve!(cache)
sol1.value ≈ 2sol2.value
```

`true`