---
# System prepended metadata

title: 'Derivatives w.r.t \epsilon_{\a'

---

### 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:

![](https://hackmd.io/_uploads/SyLdx7gza.png)

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: 
![](https://hackmd.io/_uploads/SyQwVQgG6.png)
### 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.
