This is a tiny modern Fortran library to provide utilities for Wannier interpolation. The library is meant to serve as a building block for codes that compute the resolution of quantum mechanical operators in the Brillouin zone of a crystal.

The Wannierization procedure consists of post-processing the results of a DFT calculation into Wannier states. These states are defined on real space and provide an obvious basis to express Bloch states in the reciprocal space. Many quantities in solid state physics (Berry curvature, optical responses, transport properties...) are expressed in reciprocal space in terms of the Hamiltonian and Berry connection operators and its derivatives. This library computes, given a matrix representation of the Hamiltonian and Berry connection in real space, $H_{nm}(\textbf{R})$ and $A^j_{nm}(\textbf{R})$ respectively, the matrix representation $H_{nm}(\textbf{k})$ and $A^j_{nm}(\textbf{k})$ and its derivatives in reciprocal space,

$$ H_{nm}^{lp\cdots }(\textbf{k}) = \sum_{\textbf{R}}\left[\left(iR^l\right)\left(iR^p\right)\cdots\right] e^{i\textbf{k}\cdot \textbf{R}}H_{nm}(\textbf{R}), $$

$$ A_{nm}^{j \ lp\cdots }(\textbf{k}) = \sum_{\textbf{R}}\left[\left(iR^l\right)\left(iR^p\right)\cdots\right] e^{i\textbf{k}\cdot \textbf{R}}A^j_{nm}(\textbf{R}). $$

The actual quantities $H_{nm}(\textbf{R})$ and $A^j_{nm}(\textbf{R})$ are given as inputs to the API's routines. For tight-binding models these can be given by hand. For real materials, these can be read from a file generated by the WANNIER90 code, by setting

write_tb = .true.

in the WANNIER90 file input card.

This library uses MAC's containers to store the result of the Hamiltonian and Berry connection calculations when derivatives are requested.

The derived type

type, public :: crystal

is defined.

A crystal is initialized by calling (first way)

call Cr%construct(name, &
                  direct_lattice_basis, &
                  num_bands, &
                  R_points, &
                  tunnellings, &
                  dipoles, &
                  fermi_energy)

or (second way)

call Cr%construct(name, &
                  from_file, &
                  fermi_energy)

where

• character(len=*), intent(in) :: name is the name of the crystal.
• real(dp), intent(in) :: direct_lattice_basis(3, 3) is the basis of the Bravais lattice in $\text{A}$. The first index represents the vector label and the second is the vector component such that direct_lattice_basis(i, :) is the Bravais lattice vector $\textbf{a}_i$.
• integer, intent(in) :: num_bands is the number of bands in the crystalline system.
• integer, intent(in) :: R_points(num_r_points, 3) represents the set of direct lattice points to be included in the sums $\sum_{\textbf{R}}$ in terms of the Bravais lattice basis vectors. R_points(i, :) is identified with the integers $(c^i_1, c^i_2, c^i_3)$ such that

$$ \textbf{R} = c^i_1\textbf{a}_1 + c^i_2\textbf{a}_2 + c^i_3\textbf{a}_3 $$

• complex(dp), intent(in) :: tunnellings(num_r_points, num_bands, num_bands) represents the Hamiltonian operator's matrix elements $H_{nm}(\textbf{R})$ in $\text{eV}$. The first index identifies $\textbf{R}$ and the second and third $n, m$, respectively.
• complex(dp), intent(in) :: dipoles(num_r_points, num_bands, num_bands, 3) represents the Berry connection's matrix elements $A^j_{nm}(\textbf{R})$ in $\text{A}$. The first index identifies $\textbf{R}$, the second and third $n, m$, respectively and the forth the cartesian component $j$.
• real(dp), optional, intent(in) :: fermi_energy represents the crystal's Fermi energy in $\text{eV}$.
• character(len=*), intent(in) :: from_file is the path of a WANNIER90 format tight-binding file, relative to the program's execution directory. The Bravais lattice basis, the number of bands, the number of Bravais lattice points, their coordinates and the matrix elements of the Hamiltonian and Berry connection will be read from file.

These routines retrieve crystal's properties.

Is called as,

name = Cr%name()

where character(len=120) :: name.

Is called as,

direct_lattice_basis = Cr%direct_lattice_basis()

where real(dp) :: direct_lattice_basis(3, 3) is the direct lattice basis given in the constructor in $\text{A}$.

Is called as,

reciprocal_lattice_basis = Cr%reciprocal_lattice_basis()

where real(dp) :: reciprocal_lattice_basis(3, 3) is the reciprocal lattice basis in $\text{A}^{-1}$. It obeys dot_product(reciprocal_lattice_basis(i, :), direct_lattice_basis(j, :)) $=2\pi \delta_{ij}$.

Is called as,

metric_tensor = Cr%metric_tensor()

where real(dp) :: metric_tensor(3, 3) is the metric tensor in in $\text{A}^2$. It obeys metric_tensor(3, 3) = dot_product(direct_lattice_basis(i, :), direct_lattice_basis(j, :)).

Is called as,

cell_volume = Cr%cell_volume()

where real(dp) :: cell_volume is the cell volume in $\text{A}$.

Is called as,

num_bands = Cr%num_bands()

where integer :: num_bands is the number of bands.

Is called as,

nrpts = Cr%nrpts()

where integer :: nrpts is the number of $\textbf{R}$ points.

Is called as,

rpts = Cr%rpts()

where integer :: rpts(Cr%nrpts(), 3) is the set of $\textbf{R}$ points that are included in the sums $\sum_{\textbf{R}}$ in terms of the Bravais lattice basis vectors.

Is called as,

deg_rpts = Cr%deg_rpts()

where integer :: deg_rpts(Cr%nrpts()) is the set degeneracies for each of $\textbf{R}$ points that are included in the sums $\sum_{\textbf{R}}$. For crystals constructed from input, the array is a constant 1, but for real materials constructed from file, it stores the number of sites that are equivalent.

Is called as,

hamiltonian_elements = Cr%get_real_space_hamiltonian_elements()

where complex(dp) :: hamiltonian_elements(Cr%nrpts(), Cr%num_bands(), Cr%num_bands()) is $H_{nm}(\textbf{R})$ in $\text{eV}$.

Is called as,

berry_conn = Cr%get_real_space_position_elements()

where complex(dp) :: berry_conn(Cr%nrpts(), Cr%num_bands(), Cr%num_bands(), 3) is $A^j_{nm}(\textbf{R})$ in $\text{A}$.

Is called as,

initialized = Cr%initialized()

where logical :: initialized is .true. if the crystal has been initialized and .false. otherwise.

Is called as,

fermi_energy = Cr%fermi_energy()

where real(dp) :: fermi_energy is the Fermi energy in $\text{eV}$.

Implements

$$ H_{nm}^{lp\cdots }(\textbf{k}) = \sum_{\textbf{R}}\left[\left(iR^l\right)\left(iR^p\right)\cdots\right] e^{i\textbf{k}\cdot \textbf{R}}H_{nm}(\textbf{R}), $$

in $\text{eV} \cdot (\text{A}^{d})$, where $d$ is the derivative order and is called as

H = Cr%hamiltonian(kpt [, derivative, all])

where

• real(dp), intent(in) :: kpt(3) is a point in the Brillouin zone in crystal coordinates (relative to reciprocal basis vectors), i.e., $\textbf{k}_i \in [-0.5, 0.5]$ where the Hamiltonian or its derivatives are to be computed.
• integer, optional, intent(in) :: derivative: non-negative integer. If present, is the order of the derivative of the Hamiltonian, 0 means no derivative.
• logical, optional, intent(in) :: all if present and true, all the derivatives of the Hamiltonian from 0 up to derivative are computed.
• H is the output value.
  • If derivative and all are not present, H is complex(dp) :: H(Cr%num_bands(), Cr%num_bands()) and represents the Hamiltonian $H_{nm}(\textbf{k})$.
  • If derivative = n is present and all is not present, H is MAC's complex_dp type(container) :: H and stores the Hamiltonians $n$ th derivative $H_{nm}^{lp\cdots }(\textbf{k})$. The container represents an array with shape (Cr%num_bands(), Cr%num_bands(), 3, ..., 3) where the number of indices is $n + 2$.
  • If derivative = n is present and all is present and true, H is MAC's complex_dp type(container), allocatable :: H(:) and stores, in each index, the Hamiltonians $n$ th derivative $H_{nm}^{lp\cdots }(\textbf{k})$. Each container H(i) represents an array with shape (Cr%num_bands(), Cr%num_bands(), 3, ..., 3) where the number of indices is $n + 2$. The Hamiltonian (no derivative) is stored in H(1).
  • If derivative = n is present and all is present and false, H is MAC's complex_dp type(container), allocatable :: H(:) and stores, in the first index if the container, the Hamiltonians $n$ th derivative.

Implements

$$ A_{nm}^{j \ lp\cdots }(\textbf{k}) = \sum_{\textbf{R}}\left[\left(iR^l\right)\left(iR^p\right)\cdots\right] e^{i\textbf{k}\cdot \textbf{R}}A^j_{nm}(\textbf{R}). $$

in $\text{A} \cdot (\text{A}^{d})$, where $d$ is the derivative order and is called as

A = Cr%berry_connection(kpt [, derivative, all])

where

• real(dp), intent(in) :: kpt(3) is a point in the Brillouin zone in crystal coordinates (relative to reciprocal basis vectors), i.e., $\textbf{k}_i \in [-0.5, 0.5]$ where the Berry connection or its derivatives are to be computed.
• integer, optional, intent(in) :: derivative: non-negative integer. If present, is the order of the derivative of the Berry connection, 0 means no derivative.
• logical, optional, intent(in) :: all if present and true, all the derivatives of the Berry connection from 0 up to derivative are computed.
• A is the output value.
  • If derivative and all are not present, A is complex(dp) :: H(Cr%num_bands(), Cr%num_bands(), 3) and represents the Berry connection $A_{nm}^{j}(\textbf{k})$.
  • If derivative = n is present and all is not present, A is MAC's complex_dp type(container) :: A and stores the Berry connection $n$ th derivative $A_{nm}^{j \ lp\cdots }(\textbf{k})$. The container represents an array with shape (Cr%num_bands(), Cr%num_bands(), 3, 3, ..., 3) where the number of indices is $n + 3$.
  • If derivative = n is present and all is present and true, A is MAC's complex_dp type(container), allocatable :: A(:) and stores, in each index, the Berry connection's $n$ th derivative $A_{nm}^{j \ lp\cdots }(\textbf{k})$. Each container A(i) represents an array with shape (Cr%num_bands(), Cr%num_bands(), 3, 3, ..., 3) where the number of indices is $n + 3$. The Berry connection (no derivative) is stored in A(1).
  • If derivative = n is present and all is present and false, A is MAC's complex_dp type(container), allocatable :: A(:) and stores, in the first index if the container, the Berry connection's $n$ th derivative.

The library includes an approximation of a Dirac delta as a narrow Gaussian,

$$ \delta(x) \approx \frac{1}{\sigma\sqrt{2\pi}}e^{-x^2/(2\sigma^2)}. $$

It can be called by

delta = dirac_delta(x, smr)

where

• real(dp), intent(in) :: x, smr are the evaluation point $x$ and smearing $\sigma$, respectively.
• real(dp) :: delta is an approximation to $\delta(x)$.

The library includes a utility to calculate the degeneracy of subspaces. It can be called by

dl = deg_list(eig, degen_thr)

where

• real(dp), intent(in) :: eig(n) is a list of $n$ eigenvalues $\varepsilon_i$ of an Hermitian operator sorted in ascending order.
• real(dp), intent(in) :: degen_thr is a degeneracy threshold $\lambda$ such that two eigenvalues $i, j$ obeying $|\varepsilon_i- \varepsilon_j| < \lambda$ will be considered degenerate.
• integer :: dl(n) is a list of $n$ numbers expressing the dimensionality of each subspace. If dl(i) = M, then $\varepsilon_i = \varepsilon_{i+1} = \cdots = \varepsilon_{i + M - 1}$ up to $\lambda$. If $i < j < i + M - 1$, then dl(j) = 0. If the value is nondegenerate, then dl(j) = 1.

The library includes a wrapper of dsyev and zheev routines from LAPACK. It can be called by

call diagonalize(matrix, P [, D, eig])

where

• real/complex(dp), intent(in) :: matrix(n, n) is a square symmetric/Hermitian matrix.
• real/complex(dp), intent(out) :: P(n, n) is an orthogonal/unitary matrix such that $D = P^{-1}MP$ is diagonal, $M$ being matrix.
• real/complex(dp), optional, intent(out) :: D(n, n) is the diagonal form $D$ of matrix.
• real(dp), optional, intent(out) :: eig(n) are the eigenvalues of matrix.

The library includes a utility to invert matrices by employing singular value decomposition. It can be called by

inv = inverse(matrix)

where

• real/complex(dp), intent(in) :: matrix(n, n) is an invertible square matrix $A$.
• real/complex(dp) :: inv(n, n) is an invertible square matrix $B$ obeying $AB=I$.

The library includes a wrapper of zgees routine from LAPACK. It can be called by

call schur(matrix, Z [, S, T])

where

• complex(dp), intent(in) :: matrix(n, n) is a square matrix.
• complex(dp), intent(out) :: Z(n, n) is a unitary matrix such that $S = Z^{-1}MZ$ is upper-triangular, $M$ being matrix.
• complex(dp), optional, intent(out) :: S(n, n) is the upper-triangular form $S$ of matrix.
• complex(dp), optional, intent(out) :: T(n) are the diagonal entries of S.

The library includes a wrapper of zgesvd routine from LAPACK. It can be called by

call SVD(matrix, U, V [, sigma, eig])

where

• complex(dp), intent(in) :: matrix(m, n) is a matrix.
• complex(dp), intent(out) :: U(m, m), V(n, n) are unitary matrices such that $\Sigma = U^{-1}MV$ is diagonal, $M$ being matrix.
• complex(dp), optional, intent(out) :: sigma(m, n) is the diagonal form $\Sigma$ of matrix.
• real(dp), optional, intent(out) :: eig(n) are the diagonal entries of sigma.

The library includes tools to compute the exponential of a skew-Hermitian matrix and the logarithm of a unitary matrix, expsh and logu, respectively. These can be called by

E = expsh(SH)
L = logu(U)

where

• complex(dp), intent(in) :: SH(n, n) is a skew-Hermitian ($SH^{\dagger} = -SH$) matrix.
• complex(dp) :: E(n, n) is a unitary matrix such that $E = e^{SH}$.
• complex(dp), intent(in) :: U(n, n) is a unitary matrix.
• complex(dp) :: L(n, n) is a skew-Hermitian matrix such that $L = \text{log}(U)$.

An automated build is available for Fortran Package Manager users. This is the recommended way to build and use WannInt in your projects. You can add WannInt to your project dependencies by including

[dependencies]
WannInt = { git="https://github.com/irukoa/WannInt.git" }

in the fpm.toml file.

MAC's objects

type, public :: container_specifier
type, extends(container_specifier), public :: container

are made public by WannInt.