finish writing README.

README.md

@@ -7,14 +7,14 @@ This is a tiny modern Fortran library to provide utilities for [Wannier interpol
 
 ### Working principle
 
-The [Wannierization](https://link.aps.org/doi/10.1103/RevModPhys.84.1419) procedure consists on post-processing the results of a [DFT](https://en.wikipedia.org/wiki/Density_functional_theory) calculation into Wannier states. These states are defined on real space, and provide an obvious basis to express the reciprocal basis Bloch states. 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,
+The [Wannierization](https://link.aps.org/doi/10.1103/RevModPhys.84.1419) procedure consists of post-processing the results of a [DFT](https://en.wikipedia.org/wiki/Density_functional_theory) 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_{R}\left[\left(iR^l\right)\left(iR^p\right)\cdots\right] e^{i\textbf{k}\cdot \textbf{R}}H_{nm}(\textbf{R}),
+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_{R}\left[\left(iR^l\right)\left(iR^p\right)\cdots\right] e^{i\textbf{k}\cdot \textbf{R}}A^j_{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](https://wannier.org/) code, by setting
@@ -36,6 +36,238 @@ type, public :: crystal
 ```
 is defined.
 
+## `type(crystal) :: Cr`
+### Constructor.
+
+A crystal is initialized by calling (first way)
+
+```fortran
+call Cr%construct(name, &
+                  direct_lattice_basis, &
+                  num_bands, &
+                  R_points, &
+                  tunnellings, &
+                  dipoles, &
+                  fermi_energy)
+```
+or (second way)
+```fortran
+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{\r{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{\r{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.
+
+### Handles.
+
+These routines retrieve crystal's properties.
+
+#### Name
+
+Is called as,
+
+```fortran
+name = Cr%name()
+```
+
+where `character(len=120) :: name`.
+
+#### Direct lattice basis
+
+Is called as,
+
+```fortran
+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{\r{A}}$..
+
+#### Reciprocal lattice basis
+
+Is called as,
+
+```fortran
+reciprocal_lattice_basis = Cr%reciprocal_lattice_basis()
+```
+
+where `real(dp) :: reciprocal_lattice_basis(3, 3)` is the reciprocal lattice basis in $\text{\r{A}}^{-1}$. It obeys `dot_product(reciprocal_lattice_basis(i, :), direct_lattice_basis(j, :))` $=2\pi \delta_{ij}$.
+
+#### Metric tensor
+
+Is called as,
+
+```fortran
+metric_tensor = Cr%metric_tensor()
+```
+
+where `real(dp) :: metric_tensor(3, 3)` is the metric tensor in in $\text{\r{A}}^2$. It obeys `metric_tensor(3, 3) = dot_product(direct_lattice_basis(i, :), direct_lattice_basis(j, :))`.
+
+#### Cell volume
+
+Is called as,
+
+```fortran
+cell_volume = Cr%cell_volume()
+```
+
+where `real(dp) :: cell_volume` is the cell volume in $\text{\r{A}}$.
+
+#### Number of bands
+
+Is called as,
+
+```fortran
+num_bands = Cr%num_bands()
+```
+
+where `integer :: num_bands` is the number of bands.
+
+#### Number of $\textbf{R}$ points
+
+Is called as,
+
+```fortran
+nrpts = Cr%nrpts()
+```
+
+where `integer :: nrpts` is the number of $\textbf{R}$ points.
+
+#### $\textbf{R}$ points
+
+Is called as,
+
+```fortran
+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.
+
+#### Site degeneracy
+
+Is called as,
+
+```fortran
+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.
+
+#### Hamiltonian matrix elements.
+
+Is called as,
+
+```fortran
+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}$.
+
+#### Berry connection matrix elements.
+
+Is called as,
+
+```fortran
+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{\r{A}}$.
+
+#### Fermi energy.
+
+Is called as,
+
+```fortran
+initialized = Cr%initialized()
+```
+
+where `logical :: initialized` is `.true.` if the crystal has been initialized and `.false.` otherwise.
+
+#### Initialization query.
+
+Is called as,
+
+```fortran
+fermi_energy = Cr%fermi_energy()
+```
+
+where `real(dp) :: fermi_energy` is the Fermi energy in $\text{eV}$.
+
+### Hamiltonian.
+
+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}$ and is called as
+
+```fortran
+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.
+
+### Berry connection.
+
+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{\r{A}}$ and is called as
+
+```fortran
+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.
+
+## Diagonalization utility.
+
+The library includes a wrapper of `dsyev` and `zheev` routines from LAPACK. It can be called by
+
+```fortran
+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`.
+
 # Build
 
 An automated build is available for [Fortran Package Manager](https://fpm.fortran-lang.org/) users. This is the recommended way to build and use WannInt in your projects. You can add WannInt to your project dependencies by including