### Stress definition
$$ \sigma_{\alpha\beta} = \frac{1}{\Omega}\frac{\partial E}{\partial \epsilon_{\alpha\beta}} $$

### Derivatives w.r.t $\epsilon_{\alpha,\beta}$

* volume term $\Omega$ $$ \Omega \rightarrow det(\hat{I}+\hat{\epsilon})\Omega \simeq (1+tr(\hat{\epsilon}))\Omega$$ $$ \frac{\partial\Omega}{\partial\epsilon_{\alpha,\beta}}=\delta_{\alpha\beta} \Omega$$
$$\frac{\partial{\frac{1}{\Omega}}}{\partial \epsilon_{\alpha,\beta}} = -\frac{1}{\Omega^2}\cdot \Omega \delta_{\alpha,\beta} = -\frac{\delta_{\alpha,\beta}}{\Omega}$$
* Reciprocal vectors $\mathbf{G}$ $$\mathbf{G} \rightarrow (1-\hat{\epsilon})\cdot \mathbf{G}$$ $$ \frac{\partial \mathbf{G}}{\partial\epsilon_{\alpha\beta}} = -G_{\beta}\hat{\mathbf{v}}_\alpha$$
$$ \frac{\partial G^2}{\partial\epsilon_{\alpha\beta}} = -2 G_{\alpha}G_{\beta}$$

### EXX energy with Gygi-Baldereschi correction:

Implementation briefly outlined in section B of PRB **79**,205114. Taking the notations from there the functions are: $$\phi_{\mathbf{k},\nu}(\mathbf{r}) $$; we define the codensities as:
$$ \rho^{\mathbf{k}^\prime,\nu^\prime}_{\mathbf{k},\nu}(\mathbf{r})= \phi^{~}_{\mathbf{k},\nu}(\mathbf{r})\phi^*_{\mathbf{k}^\prime,\nu^\prime}(\mathbf{r}) $$
in reciprocal space the codensities' Fourier transforms are:
$$ \rho^{\mathbf{k+q},\nu^\prime}_{\mathbf{k},\nu}(\mathbf{q+G}) \;\; \; \!\;\;\;\;\mathbf{q}=\mathbf{k-k^\prime}$$

We make first the sum on the regular MP k-point mesh and over all occupied states of the squares of the codensities
$$ A(\mathbf{q+G}) = \frac{1}{N} \sum_{\mathbf{k}}{\sum_{\nu,\nu^\prime}{\left | \rho^{\mathbf{k+q},\nu^\prime}_{\mathbf{k},\nu}(\mathbf{q+G})\right |^2}} $$

The exchange energy should then be given by the sum:
$$ \frac{4\pi e^2}{N_{\bf q}\Omega}\sum_{\bf q+G} {\frac{A({\bf q+G})}{\left | {\bf q+G}\right |^2}}$$

to correct the divergence for $\left|{\bf{q+G}}\right| = 0$ with the Gygi-Baldereschi method one has to formally add and subract to each term
$$ \frac{A(0)}{\left| \bf q+G \right|^2}e^{-\alpha\left| {\bf q+G}\right|^2}$$

This is equivalent to compute the sum above excluding ${\bf q+G}=0$ (the prime sign in the sums below) and adding the spherical averaged limit of $$\lim_{{\bf q}\rightarrow 0}{\frac{A({\bf q})-A(0)}{{\bf q}^2}} $$

plus the explicit calculation of the divergent added term:


The spherical limit instead is estimated by the difference between the sum of the regular term in dense grid that includes all ${\bf q} \neq {\bf 0}$ points minus a grid that is one half coarser in each direction. Assuming both grids are equivalent far from the center and considering that the regular grid is 8 times denser the limit can be evaluated with:
### Implementation:

#### EXX energy `exxenergy2_k`

The calculation of the energy with the formulas above is done in the routine `exxenergy2_k` in `exx.f90`. The routine computes the $A({\bf q+G})$ terms and sums them according the scheme above. The divergence term (the $D$ of formula above is associated to the $A(0)$ term.). All these term are multiplied times the *coulomb factors* and then summed to yield the total energy.

:red_circle: A factor $1/\Omega$ is added to each codensity and as we sum together the square codensities we multiply times $\Omega$ to obtain the overall $1/\Omega$ factor of the formula above for the energy.

#### Coulomb factors

In the `exx_base.f90:g2_convolution` routine computes the *coulomb factors* that in the case of HSE functional `erfc_scrlen > 0` is:
$$
\frac{4\pi e^2}{\left | {\bf q+G}\right|^2}\left( 1 - e^{-\frac{\left |{\bf q+G}\right |^2}{4 \lambda^2}}\right)
$$
and multiplying by $8/7$ (`grid_factor`) the contributions from the $\bf(q+G)$s not belonging the coarse grid and by $0$ those from the coarse grid.

#### Divergence
The divergence term includes all the terms that multiply $A(0)$ and assigned as Coulomb factor of the ${\bf q+G} = 0$ term.

These terms are calculated in the `exx_divergence` routine in `exx_base.f90`.

The routine computes explicitely the sum $$\frac{4 \pi e^2}{N_{\bf q}}\sum_{\bf q,G}^{{\bf q+G}\neq 0}\frac{fac \cdot e^{-\alpha\left|{\bf q+G}\right|^2}}{\left|{\bf q+G}\right|^2}$$

where $fac$ term can includes other terms dependent on $({\bf q+G})^2$. In the case of HSE functional $fac = 1 - e^{-q^2/4\lambda^2}$ and must be taken into account when computing the stress. This sum is done multiplying each term by `grid_factor` as above.

:red_circle: The additional terms that depend on $\alpha$, plus in other cases by the other screening parameters, are computed in a funny way with a numerical integral, they can probably be computed directly with mathematica and hard coded in the routine.

:red_circle: Other thing that need to be checked is the $\Omega$ factor added to the `aa` term.

#### Stress

Contribution to the stress for the HSE case:

* ${\bf q+G} \neq 0$ (on double grid) $$\frac{1}{\Omega^2}\left[-\frac{A({\bf q+G})}{ |{\bf q+G}|^2}\delta_{\alpha\beta} - \frac{A({\bf q+G})}{|{\bf q+G}|^4}\left( 1 - e^{-\frac{|\bf q+G|^2}{4\lambda^2}}+\frac{|{\bf q+G}|^2}{4\lambda^2} e^{-\frac{|{\bf q+G|^2}}{4\lambda^2}}\right)2(q+G)^\alpha(q+G)^\beta
\right ]
$$

* ${\bf q+G} = 0$ $$ -\frac{D A(0)}{\Omega^2}\delta_{\alpha\beta} +
\\ + \frac{A(0)}{\Omega^2} \sum_{q+G \neq 0}\frac{2(q+G)^\alpha(q+G)^\beta}{|{\bf q+G}|^4}
e^{-\alpha|{\bf q+G}|^2}{ \left[ 1 +\alpha|q+G|^2-e^{-\frac{|q+G|^2}{4\lambda^2}}\left ( 1+(\alpha+\frac{1}{4\lambda^2})|q+G|^2\right ) \right]} $$

:red_circle: the diagonal term is already implemented, we just need to check the sign and the exclusion of the volume independent term. The second term in the sum must be implemented.