Treatment Of Affine Meshes

# Motivation

Affine coordinate systems are indispensable in description of bulk systems. This is because they allow to describe a periodic system with a small unit cell containing small number of atoms [1]. Our interest to the affine coordinate system originates in the calculation of electronic density in a periodic system on an affine grid.

# Definition

Affine grid is a real space grid the points of which can be addressed with three integer numbers. Cartesian coordinate of a point in the affine grid can be defined by

(1)
\begin{align} \mathbf{r}(i,j,k) = \mathbf{u}_1 \cdot i + \mathbf{u}_2 \cdot j + \mathbf{u}_3 \cdot k + \mathbf{r}_0, \end{align}

where $\{i,j,k\}$ a triple of integers, $\{\mathbf{u}_i\}$ is a set of vectors along axes of the affine system and $\mathbf{r}_0$ is introduced to control the position of the grid. The length of the vectors determines the discretization (maximal resolution) of the grid. The equation can be written shorter by using the matrix-vector product

(2)
\begin{align} \mathbf{r}_i = \mathbf{M}_{ij} \mathbf{d}_j + \mathbf{r}_0, \end{align}

where $\mathbf{M}_{ij}$ is a matrix which is formed from the vectors $\mathbf{u}_i$ and $\mathbf{d}_j$ is the vector of integers. In the last equation and below we understand a sum over repeating indices unless the index appear on the left-hand-side of the equation.

(3)
\begin{align} \mathbf{M}_{ij} = \mathbf{u}_{ji}. \end{align}

Using the equation (2) we can find an inverse expression of the affine coordinates $\mathbf{d}_i$ in terms of the Cartesian coordinates $\mathbf{r}_i$

(4)
\begin{align} \mathbf{d}_i = \overline{\mathbf{M}}_{ij} (\mathbf{r}_j - \mathbf{r}_{0j}), \end{align}

where the overline denotes the matrix inverse.

# Gradient in affine coordinate system

(5)
\begin{align} \mathbf{\nabla}_i F(\mathbf{r}) = \frac{\partial F(\mathbf{r})}{\partial \mathbf{r}_i} \end{align}

Using the chain rule for partial derivatives and the equation (4) we can find for the Cartesian components (5)

(6)
\begin{align} \mathbf{\nabla}_i F(\mathbf{r}) = \frac{\partial F(\mathbf{r})}{\partial \mathbf{d}_j} \frac{\partial \mathbf{d}_j}{\partial \mathbf{r}_i} = \frac{\partial F(\mathbf{r})}{\partial \mathbf{d}_j} \overline{\mathbf{M}}_{ji}. \end{align}

The partial derivatives along the affine coordinates $\mathbf{d}_j$ is easy to define by any finite-element expression, and calculate the gradient in Cartesian coordinates by a vector-matrix multiplication.

# Divergence of the expressions involving gradient in affine coordinate system

## Problem posing

The problem to be solved arises in calculation of first derivatives of the exchange-correlation energy in terms of particle density of the GGA functional which can be expressed as following

(7)
\begin{align} E_{\mathrm{xc}} = \int \epsilon(n(\mathbf{r}), \mathbf{\nabla}n(\mathbf{r})) d^3 r, \end{align}

where $E_{\mathrm{xc}}$ is the exchange-correlation energy, $\epsilon(n(\mathbf{r}), \mathbf{\nabla}n(\mathbf{r}))$ is the energy density, where $n(\mathbf{r})$ is particle density and $\mathbf{\nabla}n(\mathbf{r})$ is the gradient of the particle density. Because the exchange-correlation energy depends not only on the particle density $n(\mathbf{r})$ but also on it's gradient $\mathbf{\nabla}n(\mathbf{r})$, the variational derivative of the exchange-correlation energy with respect to the particle density $v_{\mathrm{xc}}(\mathbf{r}) = \frac{\delta E_{\mathrm{xc}}}{\delta n(\mathbf{r})}$ is given by (see Euler-Lagrange equation)

(8)
\begin{align} v_{\mathrm{xc}}(\mathbf{r}) = \frac{\partial \epsilon}{\partial n} - \mathbf{\nabla} \frac{\partial \epsilon}{\partial \mathbf{\nabla} n}. \end{align}

The partial derivative of the energy density with respect to the particle density $\frac{\partial \epsilon}{\partial n}$ can be computed once a particular approximation of the exchange-correlation energy is defined. For instance a software library LIBXC delivers the values of this partial derivative for a given density and its gradient. However, all the GGA functionals indeed depend not on the vector of the gradient, but on its length and there is a tradition in quantum chemistry literature to discuss the dependency of the energy density not on the vector of the gradient of the particle density but on the length of this vector. For instance, the LIBXC library requires as input the particle density $n(\mathbf{r})$ and the square of the length of the gradient of the particle density $\sigma=|\mathbf{\nabla}n(\mathbf{r})|^2$. Moreover, the library LIBXC provides a partial derivative of the energy density with respect to the square of the gradient $\frac{\partial \epsilon}{\partial \sigma}$, rather than the partial derivative of the energy density on the vector gradient itself. In order to derive a practical expression for calculation of the exchange-correlation potential (8), we will apply a chain rule

(9)
\begin{align} v_{\mathrm{xc}}(\mathbf{r}) = \frac{\partial \epsilon}{\partial n} - 2\mathbf{\nabla} \left( \frac{\partial \epsilon}{\partial \sigma} {\mathbf{\nabla} n(\mathbf{r})}\right ). \end{align}

## Affine coordinates for GGA exchange-correlation potential

Let's derive an expression for exchange-correlation potential (9) that is ready for a finite-difference approximation when the gradient is computed on the affine grid. The function to work with we will define by the analogy to the equation (9)

(10)
\begin{align} g(\mathbf{r}) = \mathbf{\nabla}_i \left( \frac{\partial \epsilon}{\partial \sigma} {\mathbf{\nabla}_i n(\mathbf{r})}\right ). \end{align}

Let's transform the $\mathbf{\nabla}_i$ into partial derivatives along affine coordinates according to equation (6)

(11)
\begin{align} g(\mathbf{r}) = \frac{\partial}{\partial \mathbf{d}_j} \left( \frac{\partial \epsilon}{\partial \sigma} {\mathbf{\nabla}_i n(\mathbf{r})}\right ) \overline{\mathbf{M}}_{ji}. \end{align}

Thus a finite-difference expression corresponding to the last equation can be

(12)
\begin{align} g(\mathbf{r}) \approx \frac{1}{2} \left( \frac{\partial \epsilon(\mathbf{d}_j+1)}{\partial \sigma} \overline{\mathbf{M}}_{j,i} {\mathbf{\nabla}_i n(\mathbf{d}_j+1)} - \frac{\partial \epsilon(\mathbf{d}_j-1)}{\partial \sigma} \overline{\mathbf{M}}_{j,i} {\mathbf{\nabla}_i n(\mathbf{d}_j-1)}\right ). \end{align}

Probable a right thing to do for the computation of the last expression would be to organize a loop over indices $j$ in which the result will be accumulated. Moreover, in order to avoid an excessive calling of the library subroutines due to the finite-differences semi-local expression, one more quantity should be discretized on the affine grid. Namely, together with the particle density $n(\mathbf{r})$, it's gradient $\mathbf{\nabla}n(\mathbf{r})$, the values of the partial derivative $\frac{\partial \epsilon(\mathbf{r})}{\partial \sigma}$ shall be put on an affine grid. Moreover, it seems to be impossible to shrink this set of in total five quantities by using the quantum chemical conventions on GGA functionals (that energy density depends only on the lenght of the gradient of the particle density).

# Relevant modules in the source code

The set of modules is a subject of change, but one may try to look into the following modules

Modules in which the affine discretization is realized are m_dens_mesh_3d, m_mesh_3d_equ and m_bulk_vhxc_sv.

Modules in which the calculation of the GGA potential is realized are m_xc, m_numint_mol_vxc.

Bibliography
1. C. Kittel, Introduction to Solid State Physics, Wiley.