\documentstyle [12pt]{article}
%\special{ps: landscape}
%\special{ps: xy 0.5in 0.0in}
\textwidth = 25 cm
\textheight = 16 cm
\topmargin = -1.0 cm
\oddsidemargin =0.0 cm
\parindent=0pt
\def\r{{\bf r}}
\def\d{{\bf d}}
\def\k{{\bf k}}
\def\q{{\bf q}}
\def\G{{\bf G}}
\def\R{{\bf R}}

\begin{document}
{\bf Total energy and Kohn-Sham Hamiltonian of a crystal within DFT}
\par
Let us consider a crystal with $ N\rightarrow\infty$ unit cells
of volume $\Omega $, periodically repeated, with lattice vectors
$\R$. (Pseudo--)Atoms of type $\mu$ and ionic charge $Z_\mu$
are located at $\d_\mu$ in the unit cell. The system contains 
$ N \sum_\mu Z_\mu $ electrons. Its electron states are described
by $ N $ points $\k$ in the Brillouin Zone. Assuming for simplicity
a local electron-ion potential $\hat V^\mu$:
\begin{eqnarray}
E_{tot} & = & E_{kin} + E_{ion-el} + E_{Hartree} + E_{xc} + E_{ion-ion} \\
        & = & - {\hbar^2 \over 2m} \sum_{ \k,v } 
              \int \psi^*_{ \k,v } (\r)
              \nabla ^2 \psi_{ \k,v }(\r) d\r
+ \sum_{ \k,v,\mu,\R}\int \psi^*_{ \k,v } (\r)
  \hat V^\mu(\r-\d_\mu-\R) \psi_{ \k,v } (\r) d\r \nonumber \\
       &  & \mbox{}
+ {e^2 \over 2} \int {n(\r) n(\r') \over \mid \r-\r' \mid} d\r d\r'
+ \int n(\r)\epsilon_{xc}[n(\r)]d\r 
+ {e^2 \over 2} \sum'_{\mu,\nu,\R,\R'} {Z_\mu Z_\nu \over 
\mid \d_\mu+\R-\d_\nu-\R' \mid }
\end{eqnarray}
where the electron charge density $n(\r)$ is given by
\begin{equation}
n(\r) = \sum_{\k,v} {\mid \psi_{ \k,v }(\r) \mid }^2
\end{equation}
(the sum is over the lowest $ \sum_\mu Z_\mu $ occupied states for 
a semiconductor or insulator, up to the Fermi surface for a metal).
Integrals extend on all space. The primed sum appearing in the ion-ion 
term excludes terms with $\d_\mu+\R-\d_\nu-\R'=0$. 

The Kohn-Sham equation is
\begin{equation} \biggl [
- {\hbar^2 \over 2m}  \nabla ^2 
+ \sum_{\mu,\R} \hat V^\mu(\r-\d_\mu-\R) 
+ e^2 \int { n(\r') \over \mid \r-\r' \mid } d\r'
+ V_{xc}(\r)
   \biggr ] \psi_{ \k,v }(\r) 
= \epsilon_{ \k,v } \psi_{ \k,v }(\r) \end{equation}
where the exchange-correlation potential 
$ V_{xc}(\r) = (\delta E_{xc}/\delta n(\r))$.
For the LDA case only:
\begin{equation}
E_{xc}[n(\r)] = \int n(\r) \epsilon_{xc}(n(\r)) d\r, \qquad
V_{xc}(\r) = {d \over d n} \bigl(n \epsilon_{xc}(n) \bigr)_{n=n(\r)}
\end{equation}
From the Kohn-Sham equation we obtain, by summing over
occupied states:
\begin{equation}
\sum_{ \k,v } \epsilon_{ \k,v } =
- {\hbar^2 \over 2m}  \sum_{ \k,v } \int \psi^*_{ \k,v }(\r) 
   \nabla ^2 \psi_{ \k,v }(\r) d\r
+ \sum_{ \k,v,\mu,\R} \int \psi^*_{ \k,v }(\r) 
   \hat V^\mu(\r-\d_\mu-\R) \psi_{ \k,v }(\r) 
+ e^2 \int { n(\r) n(\r') \over \mid \r-\r' \mid } 
           d\r d\r'
+ \int n(\r) V_{xc}(\r) d\r 
\end{equation}
and we can give an alternate formula for the total energy of a crystal:
\begin{equation}
E_{tot} =  \sum_{ \k,v } \epsilon_{ \k,v }
- {e^2 \over 2} \int {n(\r) n(\r') \over
            \mid \r-\r' \mid} d\r d\r'
+  \int n(\r)\left(\epsilon_{xc}(\r) - V_{xc}(\r)\right) d\r 
+ {e^2 \over 2} \sum'_{\mu,\nu,\R,\R'} {Z_\mu Z_\nu \over 
\mid \d_\mu+\R-\d_\nu-\R' \mid }
\end{equation}

{\bf Plane-wave -- Pseudopotential formalism}
\par
Let us consider the \G-space representation
of the wavefunctions:
\begin{equation} |\psi_\k> = \sum_\G \Psi(\k+\G) |\k+\G>, \qquad
   \Psi ({\bf k+G}) = <{\bf k+G} \mid \psi_\k>,\qquad
   \mid {\bf k+G} > = {1 \over \sqrt V} e^{i {\bf (k+G) r}},\end{equation}
where $ V = N \Omega $ is the volume of the crystal. With these
definitions, the normalizations are:
\begin{equation}
<\k+\G|\k+\G>'=\delta_{\G,\G'},
\qquad<\psi_\k |\psi_\k>=1 \quad\mbox{if}\quad\sum_\G |\Psi(\k+\G)|^2 = 1.
\end{equation}
Let us define the Fourier trasform for a periodic function 
$ F(\r)=\sum_\R f(\r-\R) $ as:
\begin{eqnarray}
F(\G) & = & {1\over N\Omega} \int d\r F(\r) e^{-i{\bf G r}} 
= {1 \over \Omega} \int d\r f(\r) e^{-i{\bf G r}}
= <{\bf k+G}_1 \mid F(\r) \mid {\bf k+G}_2> ~,
\qquad {\bf G = G}_1-\G_2 \\
F(\r) & = &\sum_\G F(\G) e^{i{\bf G r}}.
\end{eqnarray}
We assume non local pseudopotential of general form 
$\hat V^\mu = V_\mu(r)+\sum_i V_{\mu,i}(\r,\r')$.
The total energy per unit cell in reciprocal space is:
\begin{eqnarray}
{ E_{tot} \over N } & = & {1 \over N} {\hbar^2 \over 2m} \sum_{\k,v} 
       \sum_\G {\mid \Psi_v({\bf k+G}) \mid}^2 ({\bf k+G})^2
+ \Omega \sum_\G n^*(\G) \sum_\mu S_\mu(\G) V_\mu(\G)
+ {1 \over N} \sum_{\k,v} \sum_{\mu,i} \sum_{{\bf G,G'}} 
   S_\mu({\bf G - G'}) \times \nonumber \\ & & \mbox{} \times
    \Psi_v^*({\bf k+G}) \Psi_v({\bf k+G'}) V_{\mu,i}({\bf k+G},{\bf k+G'})
+ { \Omega \over 2} \sum_\G n^*(\G) V_{Hartree}(\G)
+ \Omega\sum_\G n^*(\G)\epsilon_{xc}(\G)d\r
+ {e^2 \over 2} \sum'_{\mu,\nu,\R} 
{Z_\mu Z_\nu \over \mid \d_\mu-\d_\nu-\R\mid }
\end{eqnarray}
where $ S_\mu(\G) = \sum_\mu e^{-i{\bf G d_\mu}}$ is the structure factor, and
\begin{equation} V_{Hartree}(\G) = 4\pi e^2 { n({\bf G)} \over \G^2 }, \qquad
   V_\mu(\G)  = {1 \over \Omega} \int V_\mu(\r) 
                    e^{-i{\bf G r}} d\r, \qquad
   V_{\mu,i}(\k_1,\k_2) = {1 \over \Omega} \int
 e^{-i\k_1\r} V_{\mu,i}(\r,\r') e^{i\k_2\r'} d\r d\r' . 
\end{equation}
Note that we have assumed one atom of each kind. The generalization
is straightforward: the structure factor becomes 
$ S_\mu(\G) = \sum_{i_\mu} e^{-i{\bf G d_{i\mu}}}$
where $i_\mu$ runs over atoms of the same kind $\mu$.

Using eigenvalues sum, the total energy per unit cell is
\begin{equation}
{ E_{tot} \over N } = {1 \over N} \sum_{\k,v} \epsilon_{\k,v}
- { \Omega \over 2} \sum_\G n^*(\G) V_{Hartree}(\G)
+ \Omega \sum_\G n^*(\G) \left(\epsilon_{xc}(\G)-V_{xc}(\G)\right)
+ {e^2 \over 2} \sum'_{\mu,\nu,\R} 
  {Z_\mu Z_\nu \over \mid \d_\mu-\d_\nu-\R \mid }.
\end{equation}
In the plane-wave representation the Kohn-Sham equation becomes
\begin{equation}
 \sum_{\G'} <{\bf k+G} \mid H-\epsilon \mid {\bf k+G'}> \Psi ({\bf k+G'}) = 0,
\qquad\mbox{or} \qquad
 \sum_{\G'} <{\bf k+G} \mid H \mid {\bf k+G'}> \Psi ({\bf k+G'}) =
 \epsilon \Psi ({\bf k+G})
\end{equation}
The matrix elements of the hamiltonian are
\begin{eqnarray} <{\bf k+G} \mid H-\epsilon \mid {\bf k+G'}> & = &
 \biggl ( - {\hbar^2 \over 2m} ({\bf k+G})^2 -\epsilon \biggr ) \delta_{\bf GG'}
+ \sum_\mu S_\mu(\G-\G') \biggl (  V_\mu(\G-\G')
+ \sum_i V_{\mu,i}({\bf k+G},{\bf k+G'}) \biggr ) \nonumber \\
& + & V_{Hartree}(\G-\G') + V_{xc}(\G-\G') .
\end{eqnarray}

{\em Divergent Terms in the potential}
\par
The Hartree term, $V_{Hartree}(0)$, and local potential term,
$\sum_\mu S_\mu(0)V_\mu(0)$, are separately divergent and must 
be treated in a special way.
Let us consider their sum $\widetilde V(\r)=V_{loc}(\r)+V_{Hartree}(\r)$.
Its $\G=0$ term is not divergent:
\begin{equation}
    \widetilde V(\G=0) ={1\over\Omega} \int d\r\left(\sum_\mu V_\mu(\r-\d_\mu) 
        + {1\over N} e^2 \int {n(\r')\over |\r-\r'|} d\r'\right)=
       {1\over\Omega} \sum_\mu \int d\r \left( V_\mu(r) 
       + {Z_\mu e^2\over r} \right ) = {1\over\Omega} \sum_\mu \alpha_\mu
\end{equation}
where we used
\begin{equation}
    V_\mu(r) \sim -{Z_\mu e^2\over r}\quad\mbox{for large~} r,\qquad
    {1\over N} \int n(\r) = \sum_\mu Z_\mu.
\end{equation}

The $\alpha_\mu$ are parameters depending only on the pseudopotential.

\newpage
{\em Divergent Terms in the energy}

\par
The $\G=0$ terms of the ion-ion, Hartree, and local
pseudopotential terms in the total energy are separately divergent
and must be treated in a special way. Let us call $E_{div}$ the
sum of all divergent terms.

First Step: split $E_{div} = E_{div}^{(1)} + E_{div}^{(2)}$, with
\begin{equation}
E_{div}^{(1)} = \int n(\r) \sum_\mu V_\mu(\r-\d_\mu) d\r
        + {1\over N} e^2 \int {n(\r)n(\r')\over |\r-\r'|} d\r d\r'
\end{equation}
\begin{equation}
E_{div}^{(2)} = {e^2\over 2} 
      \sum'_{\mu,\nu,\R}{Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
    - {1\over N} {e^2\over 2} \int {n(\r)n(\r')\over |\r-\r'|} d\r d\r'
\end{equation}
Using the previous definition of $\widetilde V(\r)$, 
the first divergent term can be written as
\begin{equation}
E_{div}^{(1)}= \int n(\r) \widetilde V(\r) d\r.
\end{equation}
The $\G=0$ term of $\widetilde V(\G)$ is not divergent and has been
previously calculated:
\begin{equation}
 \widetilde V(\G=0) = {1\over\Omega} \sum_\mu \alpha_\mu,\qquad
 n(\G=0)=\sum_\mu {Z_\mu \over \Omega}.
\end{equation}
We finally get for the $\G=0$ contribution what is usually called
``$\alpha Z$ term'':
\begin{equation}
E_{div}^{(1)} = \Omega\sum_{\G\neq 0} n^*(\G) \widetilde V(\G)
              + {1\over\Omega}(\sum_\mu Z_\mu)(\sum_\mu \alpha_\mu)
              = \Omega\sum_\G n^*(\G) \widetilde V_{loc}(\G)
              + 2 \widetilde E_{Hartree}
\end{equation}
where $\widetilde V_{loc}$ is the local potential for $\G\neq 0$,
contains the $\alpha Z$ term in $\G=0$ component, and
\begin{equation}
\widetilde E_{Hartree} = {\Omega\over 2}\sum_{\G\neq 0} n^*(\G)V_{Hartree}(\G).
\end{equation}
Second step: write
$E_{div}^{(2)}= E_{Ewald}^{(1)} + E_{Ewald}^{(2)}-E_{Hartree}$,
with
\begin{equation}
E_{Ewald}^{(1)} = 
  {e^2\over 2} \sum'_{\mu,\nu,\R}{Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
  \mbox{erfc}(\sqrt{\eta}|\d_\mu-\d_\nu-\R|)~,\qquad
E_{Ewald}^{(2)} = 
   {e^2\over 2} \sum_{\mu,\nu,\R}{Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
   \mbox{erf} (\sqrt{\eta}|\d_\mu-\d_\nu-\R|)
   - e^2 \sqrt{\eta\over\pi}\sum_\mu Z_\mu^2.
\end{equation}
This identity is verified for any value of $\eta$.
The sum in $E_{Ewald}^{(2)}$ includes the term with $\d_\mu-\d_\nu-\R=0$
(note the missing prime), that is subtracted back in the second term of 
$E_{Ewald}^{(2)}$ (note that 
$\mbox{erf}(x)\rightarrow 2x/\sqrt\pi$ for small $x$).

The first Ewald term $E_{Ewald}^{(1)}$ is rapidly convergent 
in real space for any reasonable values of $\eta$.

The sum in $E_{Ewald}^{(2)}$ can be written
as the interaction energy between point charges $n_c(\r)$ 
and the potential $V_g(\r)$ produced by a gaussian distribution
of charges:
\begin{equation}
E_{Ewald}^{(2)}= {1\over 2}\int n_c(\r)V_g(\r)d\r
   - e^2 \sqrt{\eta\over\pi}\sum_\mu Z_\mu^2~,\qquad
n_c(\r)=\sum_\mu Z_\mu\delta(\r-\d_\mu)~,
\qquad V_g(\r)=e^2\sum_{\mu,\R} 
{Z_\mu\mbox{erf}(\sqrt{\eta}|\r-\d_\mu-\R|) \over |\r-\d_\mu-\R|}
\end{equation}
In reciprocal space, by using the Fourier transform
\begin{equation}
   {1\over r'} \mbox{erf} (\sqrt{\eta}r')
= \left({\eta\over \pi}\right)^{3/2} \int {e^{-\eta r^2}\over |\r-\r'|} d\r
 = \int {4\pi e^{-G^2/4\eta}\over G^2}e^{i\G\cdot\r'} d\G
\end{equation}
one obtains (forgetting for the moment the divergence of 
$V_g(\G=0)$):
\begin{equation}
E_{Ewald}^{(2)}= {\Omega\over 2}\sum_{\G} n_c^*(\G)V_g(\G)
   - e^2 \sqrt{\eta\over\pi}\sum_\mu Z_\mu^2~,\qquad
n_c(\G)={1\over\Omega}\sum_\mu Z_\mu e^{i\G\cdot\d_\mu}~,
\qquad V_g(\G) =  {4\pi e^2\over\Omega}
\sum_\mu Z_\mu e^{i\G\cdot\d_\mu} {e^{-G^2/4\eta}\over G^2}
\end{equation}
The $\G = 0$ contribution to $E_{Ewald}^{(2)}-E_{Hartree}$:
\begin{equation}
E_0 = {\Omega\over 2} \left(n_c(0)V_g(0)-n(0)V_{Hartree}(0)\right)
\end{equation}
is no longer divergent, because $n(0)=n_c(0)=\sum_\mu Z_\mu/\Omega$ 
due to the neutrality of the system:
\begin{eqnarray}
(V_g-V_{Hartree})(\G=0) & = & {e^2\over N \Omega}
 \int\left( \sum_{\mu,\R} Z_\mu
{\mbox{erf}(\sqrt{\eta}|\r-\d_\mu-\R)|)\over |\r-\d_\mu-\R)|}
- \int {n(\r')\over |\r-\r'|}d\r'\right)d\r \nonumber \\
& = &  {e^2\over\Omega}\left(\sum_\mu Z_\mu\right)
  \int {\mbox{erf}(\sqrt{\eta}r)-1\over r}d\r 
  = {e^2\over \Omega} \left(\sum_\mu Z_\mu\right) {\pi\over\eta}
\end{eqnarray}
The integral appearing in the last expression can be found
in tables:
\begin{equation}
  \int {\mbox{erf}(\sqrt{\eta}r)-1\over r}d\r = 
  4\pi \int (\mbox{erf}(\sqrt{\eta}r)-1)rdr = 4\pi{1 \over 4\eta}.
\end{equation}
Putting all pieces together, one obtains for $E^{(2)}_{div}$:
\begin{eqnarray}
E^{(2)}_{div} = -\widetilde E_{hartree} + E_{Ewald}
  & = & - {\Omega\over 2} \sum_{\G\neq0} n^*(\G)V_{Hartree}(\G)
                + {4\pi\over\Omega} {e^2\over 2}\sum_{\G\neq 0}
                  \left |\sum_\mu Z_\mu e^{i\G\d_\mu}\right |^2
                  {e^{-G^2/4\eta}\over G^2} \nonumber \\ & & \mbox{}
                + {e^2\over 2} \sum'_{\mu,\nu,\R}
                  {Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
                  \mbox{erfc}(\sqrt{\eta}|\d_\mu-\d_\nu-\R|)
                - e^2 \sqrt{\eta\over\pi}\sum_\mu Z_\mu^2
                - {4\pi\over\Omega} {e^2\over 2} {1\over 4\eta}
                  \left(\sum_\mu Z_\mu\right)^2
\end{eqnarray}
and for the total energy:
\begin{eqnarray}
{ E_{tot} \over N } & = & {1 \over N} {\hbar^2 \over 2m} \sum_{\k,v} 
       \sum_\G {\mid \Psi_v({\bf k+G}) \mid}^2 ({\bf k+G})^2
+ {1 \over N} \sum_{\k,v} \sum_{\mu,i} \sum_{{\bf G,G'}} 
   S_\mu({\bf G - G'}) \Psi_v^*({\bf k+G}) \Psi_v({\bf k+G'}) 
   V_{\mu,i}({\bf k+G},{\bf k+G'}) \nonumber \\ & & \mbox{} 
+ \Omega \sum_\G n^*(\G) \epsilon_{xc}(\G)
+ \Omega \sum_\G n^*(\G) \widetilde V_{loc}(\G)
+ \widetilde E_{Hartree} + E_{Ewald}.
\end{eqnarray}
One can use the sum of the eigenvalues to calculate the total energy: 
the kinetic, non local, and local (including the $\alpha Z$ term) 
contributions disappear and the expression of the total energy
becomes: 
\begin{eqnarray}
{ E_{tot} \over N } = {1 \over N} \sum_{\k,v} \epsilon_{\k,v}
               + \Omega \sum_\G n^*(\G)\left(\epsilon_{xc}(\G)-V_{xc}(\G)\right)
               - \tilde E_{Hartree} + E_{Ewald}.
\end{eqnarray}

{\em Calculation of energy in CP}

\par
Let us define $n_g$ as the sum of gaussian charges centered at 
atomic sites:
\begin{equation}
n_g(\r) = \sum_\mu Z_\mu \left ({\eta\over\pi}\right)^{3/2}
                         e^{-\eta(\r-\d_\mu)^2}.
\end{equation}
The divergent terms $E_{div}$ of the energy can be rewritten as
$E_{div} = E_{div}^{(1)} + E_{div}^{(2)}$, with
\begin{equation}
E_{div}^{(1)} = \int n(\r) \sum_\mu V_\mu(\r-\d_\mu) d\r
        - {1\over N} e^2 \int {n(\r)n_g(\r')\over |\r-\r'|} d\r d\r'
\end{equation}
\begin{equation}
E_{div}^{(2)} = {e^2\over 2} 
      \sum'_{\mu,\nu,\R}{Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
    + {1\over N} {e^2\over 2} \int {n(\r)n(\r')\over |\r-\r'|} d\r d\r'
    + {1\over N} e^2 \int {n  (\r)n_g(\r')\over |\r-\r'|} d\r d\r'
\end{equation}
$E_{div}^{(1)}$ can be rewritten as
\begin{equation}
E_{div}^{(1)}= \int n(\r) V_{loc+g}(\r)d\r,\qquad
 V_{loc+g}(\r) = \left ( \sum_\mu V_\mu(\r-\d_\mu)  + V_g(\r) \right ),
\end{equation}
where $V_g(\r)$ is the potential generated by the gaussians.
The singularity at $\G=0$ disappears :
\begin{equation}
 V_{loc+g}(\G=0) = {1\over\Omega} \sum_\mu \int d\r 
    \left(V_\mu(r) + {Z_\mu e^2\mbox{erf}(\sqrt{\eta}r)\over r}\right)
   \equiv  {1\over\Omega} \sum_\mu \alpha'_\mu.
\end{equation}
Note that this term is similar to, but not equal to, the $\alpha Z$ term.
$E_{div}^{(1)}$ can be rewritten as:
\begin{equation}
E_{div}^{(1)} =  {1\over\Omega} \left(\sum_\mu Z_\mu\right)
                                \left(\sum_\mu \alpha'_\mu\right)
        + \Omega \sum_{\G\neq 0} n^*(\G)V_{loc+g}(\G) .
\end{equation}

$E_{div}^{(2)}$ can be rewritten as $E_{div}^{(2)}= E_{Ewald} + E_{H+g}$,
where $E_{Ewald}$ is the well-known Ewald sum,
\begin{equation}
E_{Ewald} = {e^2\over 2} 
      \sum'_{\mu,\nu,\R}{Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
    - {1\over N} {e^2\over 2} \int {n_g(\r)n_g(\r')\over |\r-\r'|} d\r d\r',
\end{equation}
while $E_{H+g}$ is the electrostatic energy of a system of electrons
and ions with a gaussian charge distribution:
\begin{equation}
E_{H+g} = {1\over N} {e^2\over 2} 
         \int {(n(\r)+n_g(\r))(n(\r')+n_g(\r'))\over |\r-\r'|} d\r d\r'.
\end{equation}
Both terms are regular at $\G=0$ because they are the electrostatic energy 
of neutral systems. $E_{H+g}$ can be directly calculated:
\begin{equation}
E_{H+g} = \Omega \sum_{\G\neq 0} (n^*(\G)+n_g^*(\G)) V_{H+g}(\G), \quad
V_{H+g}(\G) = {4\pi e^2\over G^2} \left(n(\G)+n_g(\G) \right) .
\end{equation}

$E_{Ewald}$ can be easily computed using the same technique used before:
\begin{equation}
E_{Ewald} = E^{(1)}_{Ewald} + E^{(2)}_{Ewald}, 
\end{equation}
where (note the sum over all vectors including $\d_\mu-\d_\nu-\R=0$):
\begin{equation}
E^{(1)}_{Ewald} = 
  {e^2\over 2} \sum_{\mu,\nu,\R}{Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
  \mbox{erf}(\sqrt{\zeta}|\d_\mu-\d_\nu-\R|)
    - {1\over N} {e^2\over 2} \int {n_g(\r)n_g(\r')\over |\r-\r'|} d\r d\r'
\end{equation}
and (note the self-interaction term compensating the $\d_\mu-\d_\nu-\R=0$
contribution of the former term):
\begin{equation}
E^{(2)}_{Ewald} = 
  {e^2\over 2} \sum'_{\mu,\nu,\R}{Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
  \mbox{erfc}(\sqrt{\zeta}|\d_\mu-\d_\nu-\R|)
   - e^2 \sqrt{\zeta\over\pi}\sum_\mu Z_\mu^2.
\end{equation}
for an arbitrary value of $\zeta$.
$E^{(1)}_{Ewald}$ can be made to vanish, because both terms appearing 
in it can be written in reciprocal space as an Ewald sum. Leaving apart
the $\G=0$ contribution,
\begin{equation}
  {1\over N} {e^2\over 2} \int {n_g(\r)n_g(\r')\over |\r-\r'|} d\r d\r'
 = {\Omega \over 2}\sum_{\G\neq 0} { n_g^*(\G)V_g(\G)}
 = {4\pi\over\Omega}{e^2\over 2}\sum_{\G\neq 0}
    \left|\sum_\mu Z_\mu e^{i\G\cdot\d_\mu}\right|^2 {e^{-G^2/2\eta}\over G^2}.
\end{equation}
This expression comes from the Fourier transform of $n_g(\r)$ and $V_g(\r)$:
\begin{equation}
n_g(\G) = {1\over \Omega} 
\sum_\mu Z_\mu e^{i\G\cdot\d_\mu} {e^{-G^2/4\eta}},\quad
V_g(\G) = {4\pi e^2 n_g(\G) \over G^2}.
\end{equation}
By setting $2\zeta=\eta$, one finds exactly the reciprocal space
expression for the first term:
\begin{equation}
  {e^2\over 2} \sum_{\mu,\nu,\R}{Z_\mu Z_\nu\over|\d_\mu-\d_\nu-\R|}
  \mbox{erf}(\sqrt{\zeta}|\d_\mu-\d_\nu-\R|)
 = {4\pi\over\Omega}{e^2\over 2}\sum_{\G\neq 0}
    \left|\sum_\mu Z_\mu e^{i\G\cdot\d_\mu}\right|^2 {e^{-G^2/4\zeta}\over G^2}.
\end{equation}
$E^{(2)}_{Ewald}$ is a rapidly convergent sum in real space,
plus a term coming from the self-interaction of gaussians.

{\em Miscellaneous}

\par
When a set of special points $\{ \k_i \} $, with weights 
$ w_i$, $ \sum_i w_i = 1$, 
is used to sample the Brillouin Zone, one has:
\begin{equation}
{1 \over N} \sum_\k f(\k) \Longrightarrow 
\sum_i w_i f(\k_i).
\end{equation}

The $\psi(\r)$ as defined above are vanishingly small in order 
to be normalized. What is actually calculated, and used in the 
Fast Fourier Transform algorithm, is $\sqrt{N} \psi(\r)$:
$\Psi(\k+\G) \stackrel{FFT}{\longleftrightarrow}\sqrt{N}\psi(\r)$.
This ensures the correct normalization of the charge density.

\newpage
\vskip 0.5 truecm
{\bf Matrix elements of non-local pseudopotentials}
\vskip 0.5 truecm
{\em Semilocal form}:
\begin{equation} 
\hat V^\mu = V_\mu(r) + \sum_l V_{\mu,l}(r) \hat P_l
\end{equation}
where $ \hat P_l= |l><l| $ is the projector on angular momentum l.
Using the expansion of plane waves into spherical harmonics
$Y_{lm}$ and spherical Bessel functions $j_l$:
\begin{equation}
e^{-i\k\r}=(2l+1)\sum_l i^lj_l(kr)P_l(\k\cdot \r)
           = 4\pi \sum_l i^lj_l(kr)\sum_m Y_{lm}^*(\k)Y_{lm}(\r)
\end{equation}
one gets
\begin{equation} V_{\mu,l}(\k_1,\k_2) = {1 \over \Omega} \int
e^{-i \k_1 \cdot \r} V_{\mu,l}(\r)  \hat P_l e^{ i \k_2 \cdot \r} d\r
= {4\pi \over \Omega} (2l+1) P_l(\k_1\cdot \k_2)
\int_0^\infty r^2 j_l(k_1 r)j_l(k_2 r) V_{\mu,l}(r) dr
={4\pi \over \Omega} (2l+1) P_l(\k_1\cdot \k_2) F_l^\mu(k_1,k_2)
\end{equation}
where the $P_l$ are Lagrange polynomials.
Using Bachelet-Hamann-Schl\"uter parameterization:
\begin{equation}
 \hat V^\mu(\r) = -{Z_\mu e^2 \over r}\sum_{n=1}^{n_c} c_n^\mu
              \mbox{erf}(\sqrt{\alpha_n^\mu}r)
 + \sum_{l=0}^{\hat l} \sum_{n=1}^{n_l} \biggl(A_{n,l}^\mu+A_{n+3,l}^\mu r^2
   \biggr)   e^{-\alpha_{n,l}^\mu r^2}   \hat P_l
\end{equation}
where $n_c=2$,$n_l=3$, the $F_l^\mu(k_1,k_2)$ have an analytical expression.
For the $\alpha_\mu$ one finds
\begin{equation}
\alpha_\mu = \int d\r\left(V_{loc}(r)+{Z_\mu e^2\over r}\right) =
4\pi Z_\mu e^2 \int r^2 dr \left({-c_1 \mbox{erf}(\sqrt{\alpha_1}r)
-(1-c_1)\mbox{erf}(\sqrt{\alpha_2}r)\over r} -{1\over r}\right) = 
4\pi Z_\mu e^2 \left({c_1\over 4\alpha_1}+{1-c_1\over 4\alpha_2}\right).
\end{equation}

{\em Separable form}:
\begin{equation} 
\hat V^\mu = V_\mu(r) + \sum_l {1\over v_l} f_{\mu,l}(r)f^*_{\mu,l}(r').
\end{equation}
For simplicity, we consider just one term per value of $l$.
Using the same expansion as above one finds
\begin{equation} 
V_{\mu,l}(\k_1,\k_2) = {1 \over \Omega}{1\over v_l} 
\int e^{-i \k_1 \cdot \r} f_{\mu,l}(\r) d\r
\int e^{ i \k_2 \cdot \r'}f^*_{\mu,l}(\r')d\r'
= {4\pi\over \Omega} {1\over v_l}
\sum_{m=-l}^l\int r^2 j_l(k_1r)f_{\mu,l}(\r)(r)dr Y_{lm}(\hat \k_1)
             \int r^2 j_l(k_2r)f^*_{\mu,l}(\r)(r)dr Y^*_{lm}(\hat\k_2).
\end{equation}
The angular term could be summed to yield Lagrange polynomials, but
it is much more convenient to keep the ``separated'' form in
Fourier space for computational purposes.

{\bf Beyond LDA}

For spin-polarised systems, one introduces two distinct charge densities
for spin up and down, $n_\sigma$, where $\sigma=+,-$, with 
$n(\r)$ the total charge density, $n(\r)=n_+(\r)+n_-(\r)$. The 
exchange-correlation energy $E_{xc}$ is written as a functional of
both $n_\sigma$'s:
\begin{equation}
E_{xc}= \int n(\r) \epsilon_{xc}(n_+(\r),n_-(\r)) d\r.
\end{equation}
The corresponding exchange-correlation potentials $ V^\sigma_{xc}(\r)$
are given by
\begin{equation}
V^\sigma_{xc}(\r) = {\delta E_{xc}\over\delta n_\sigma(\r) }
                  = \epsilon_{xc}(n_+(\r),n_-(\r)) + n(\r) 
                    {\delta \epsilon_{xc}(n_+(\r),n_-(\r))
                      \over \delta n_\sigma(\r) }
      = \epsilon_{xc}(n_+(\r),n_-(\r))
      + n(\r)\left .{\partial \epsilon_{xc}(n_+,n_-)
                    \over \partial n_\sigma }\right|_{n_+=n_+(\r),n_-=n_-(\r)}
\end{equation}
(the functional derivatives can be performed as normal derivatives).
Approximate forms of the exchange-correlation energy are traditionally
expressed in terms of the density parameter $r_s=(3/4\pi n)^{1/3}$ 
and of the spin polarisation $\zeta=(n_+-n_-)/n$. The following 
derivatives are useful in calculating the potential:
\begin{equation}
{\partial r_s \over \partial n_\sigma} = - {r_s\over 3n}~.\qquad
{\partial\zeta\over \partial n_\sigma} = -{\zeta-\mbox{sgn}\sigma\over n}
\end{equation}
where $\mbox{sgn}\sigma=+1$ if $\sigma=+$, $\mbox{sgn}\sigma=-1$
if $\sigma=-$. 
If $\epsilon_{xc}=\epsilon_{xc}(r_s,\zeta)$, the following expression
holds for $V^\sigma_{xc}$:
\begin{equation}
V^\sigma_{xc}(r_s,\zeta) = \epsilon_{xc}(r_s,\zeta)
                         - {r_s\over 3}{\partial \epsilon_{xc}(r_s,\zeta)
                                        \over\partial r_s} 
                         - (\zeta -\mbox{sgn}\sigma)
                           {\partial \epsilon_{xc}(r_s,\zeta)
                                        \over\partial \zeta}
\end{equation}
The exchange part, $\epsilon_x$, of $\epsilon_{xc}$ 
is usually written as
\begin{equation}
\epsilon_x(r_s,\zeta)=-{3\over 4\pi r_s} \left({9\pi\over 4}\right)^{1/3}
                         {(1+\zeta)^{4/3}+(1-\zeta)^{4/3}\over 2}.
\end{equation}
The corresponding exchange potential is given by
\begin{equation}
V^\sigma_x(r_s,\zeta) = {4\over 3} \epsilon_x(r_s,\zeta)
                      + (\zeta -\mbox{sgn}\sigma)
                       {1\over \pi r_s} \left({9\pi\over 4}\right)^{1/3}
                         {(1+\zeta)^{1/3}-(1-\zeta)^{1/3}\over 2}.
\end{equation}
Note that with the above form of $\epsilon_{xc}$ the exchange energy
$E_x$ is simply
\begin{equation}
E_x[n_+,n_-]= -{3\over 4} \left({6\over \pi}\right)^{1/3} \int
[n_+(\r)^{4/3}+n_-(\r)^{4/3}] d\r.
\end{equation}

{\bf Gradient-corrected functionals}

The exchange-correlation energy $E_{xc}$ is written as a functional of
both $n_\sigma$'s and of the gradients $\nabla n_\sigma$'s:
\begin{equation}
E_{xc}=\int n(\r)\epsilon_{xc}(n_+(\r),n_-(\r),\nabla n_+,\nabla n_-)d\r.
\end{equation}
The corresponding exchange-correlation potentials $ V^\sigma_{xc}(\r)$
are given by
\begin{equation}
V^\sigma_{xc}(\r) = {\delta E_{xc}\over\delta n_\sigma(\r) }
                  = {\partial F_{xc}\over\partial n_\sigma(\r)}
                  -\sum_\alpha {\partial \over \partial \r_\alpha}
                   \left [{\partial F_{xc}\over\partial
                           (\nabla_\alpha n_\sigma) }\right]
                  = {\partial F_{xc}\over\partial n_\sigma(\r)}
                  -\sum_\alpha {\partial \over \partial\r_\alpha} 
                   \left [
                     {\partial F_{xc}\over\partial |\nabla n_\sigma|}
                     {\nabla_\alpha n_\sigma\over|\nabla n_\sigma| }
                   \right]
\end{equation}
where $F_{xc}=n\epsilon_{xc}$ and $\alpha$ denotes cartesian components.
In plane-wave codes, the derivative in the second term is usually done
using FFT.

{\em Exchange}: The gradient-corrected  exchange energy $E_x$ is 
usually written as
\begin{equation}
E_x[n_+,n_-,\nabla n_+,\nabla n_-] = 
  {1\over2} E_x[2n_+]+{1\over2} E_x[2n_-]
\end{equation}
where
\begin{equation}
E_x[n] =  \int n\epsilon_x(r_s,s)d\r=
          \int n\epsilon^0_x(r_s) F(s)d\r~,\qquad
\epsilon^0_x(r_s) = -{3k_F\over 4\pi}
\end{equation}
where $k_F=(3\pi^2n)^{1/3}=1.9191583/r_s$, $s=|\nabla n|/2k_Fn$,
$\epsilon^0_x(r_s)$ is the usual Slater exchange
and $F(s)$ is an ``enhancement factor''

Each component (up and down) of the exchange potential 
contains a first term:
\begin{equation}
{\partial\over\partial n_\sigma}(n\epsilon_x(r_s,s)) = 
\epsilon_x + n{\partial \epsilon_x\over\partial n_\sigma}+
      n{\partial \epsilon_x\over\partial s}
       {\partial s\over\partial n_\sigma}=
     {4\over 3}\epsilon^0_xF(s) + 
      n\epsilon^0_x{d F\over d s}
                   {\partial s\over\partial n_\sigma}=
     {4\over 3}\epsilon^0_x\left(F(s) - s{d F\over d s}\right).
\end{equation}
The last result uses the following identity:
\begin{equation}
 {\partial s\over\partial n_\sigma} = 
 -{|\nabla n|\over 2n^2k_F}{\partial n\over\partial n_\sigma}
 -{|\nabla n|\over 2nk_F^2}{\partial k_F\over\partial n_\sigma}
=-{|\nabla n|\over 2n^2k_F}-{|\nabla n|\over 2nk_F} {1\over 3n}
= -{4s\over 3n}.
\end{equation}
Second term of the exchange potential:
\begin{equation}
{\partial \over\partial|\nabla n_\sigma|}(n\epsilon_x(r_s,s)) =
 n {\partial \epsilon_x\over\partial s}
   {\partial s\over\partial|\nabla n_\sigma|} =
  {1\over 2k_F} {\partial\epsilon_x \over\partial s} =
  {1\over 2k_F} \epsilon^0_x  {d F\over d s} =
  -{3\over 8\pi} {d F\over d s} 
\end{equation}

{\em Correlation}: Gradient-correction contributions to the 
correlation energy $E_c$ are usually expressed as a functions
of variables $t,r_s,\zeta$:
\begin{equation}
E_c[n_+,n_-,\nabla n_+,\nabla n_-] = 
   \int n(\r) \left [\epsilon_c(r_s,\zeta)+H(t,r_s,\zeta)\right] d\r,
\end{equation}
where $\epsilon_c(r_s,\zeta)$ is the LSDA correlation
energy; $t={|\nabla n|/(2gn k_s)}$,
$g=[(1+\zeta)^{2/3}+(1-\zeta)^{2/3}]/2$, 
$k_s=(4k_F/\pi)^{1/2}$.
Note that
$k^2_s/k^2_F=0.663436444 r_s$.

Gradient-correction contribution to the correlation potential,
first term:
\begin{equation}
{\partial\over\partial n_\sigma}(nH(t,r_s,\zeta)) = 
H + n{\partial H\over\partial n_\sigma}=
H + n{\partial H\over\partial t}
     {\partial t\over\partial n_\sigma}
  + n{\partial H\over\partial r_s}
     {\partial r_s\over\partial n_\sigma}
  + n{\partial H\over\partial \zeta}
     {\partial\zeta\over\partial n_\sigma}.
\end{equation}
Note that
\begin{equation}
 {\partial t\over\partial n_\sigma} = 
 -{|\nabla n|\over 2g^2nk_s}{\partial g\over\partial n_\sigma}
 -{|\nabla n|\over 2gn^2k_s}{\partial n\over\partial n_\sigma}
 -{|\nabla n|\over 2gnk_s^2}{\partial k_s\over\partial n_\sigma}
=-{|\nabla n|\over 2gnk_s} 
  \left({\partial g\over\partial\zeta}{\partial\zeta\over\partial n_\sigma}
       +{1\over n} + {1\over 6n}\right)
= t \left({(1+\zeta)^{-1/3}-(1-\zeta)^{-1/3}\over 3}
         {\zeta-\mbox{sgn}\sigma\over n} + {7\over 6n}\right).
\end{equation}
Using preceding results for $\partial r_s/\partial n_\sigma$ and
$\partial\zeta/\partial n_\sigma$, one finds
\begin{equation}
{\partial\over\partial n_\sigma}(nH(t,r_s,\zeta)) = 
H + t\left({(1+\zeta)^{-1/3}-(1-\zeta)^{-1/3}\over 3}
         (\zeta-\mbox{sgn}\sigma) + {7\over 6}\right)
    {\partial H\over\partial t} 
  - {r_s\over 3} {\partial H\over\partial r_s}
  - (\zeta-\mbox{sgn}\sigma) {\partial H\over\partial \zeta}.
\end{equation}

Second term for the correlation potential:
\begin{equation}
{\partial \over\partial|\nabla n_\sigma|}(nH(t,r_s,\zeta)) =
 n {\partial H\over\partial t}
   {\partial t\over\partial|\nabla n_\sigma|} =
  {1\over 2gk_s} {\partial H\over\partial t}
\end{equation}


{\bf PW91 Gradient correction}

{Exchange term}: the enhancement function $F$ is given by
\begin{equation}
F(s) = {1+0.19645 s \mbox{ sinh}^{-1}(7.7956s)+(0.2743-0.1508 e^{-100s^2})s^2
  \over 1+0.19645 s \mbox{ sinh}^{-1}(7.7956s) +0.0004s^4}.
\end{equation}

{Correlation term}: the function $H$ is the sum of two terms,
$H = H_0+H_1$:
\begin{equation}
H_0 = g^3{\beta^2\over2\alpha}\mbox{log}\left(1+{2\alpha\over\beta} 
t^2 {1+At^2\over 1+At^2+A^2t^4 }\right),
\qquad
H_1 = \nu \left(C_c(r_s)-C_c(0)-{3\over 7} C_x\right) 
      g^3 t^2 e^{-a_4 (k^2_s/k^2_F) g^4 t^2}
\end{equation}
where
\begin{equation}
A = {2\alpha\over\beta} 
    {1\over e^{-2\alpha\epsilon_c(r_s,\zeta)/(g^3\beta^2)}-1}~,\quad
C_c(r_s) = -C_x + {c_1+c_2 r_s+c_3 r_s^2\over1+c_4 r_s+c_5 r_s^2+c_6 r_s^3}.
\end{equation}
The numerical factors are:
$\nu=(16/\pi)(3\pi^2)^{1/3}=15.75592$, $C_c(0)=0.004235$, 
$C_x=-0.001667212$, $\alpha=0.09$,
$c_1=0.002568$, $c_2=0.023266$, $c_3=7.389\times 10^{-6}$,
$c_4=8.723$, $c_5=0.472$, $c_6=7.389\times 10^{-2}$, 
$a_4=100$, $\beta=\nu C_c(0)$.

Let us introduce the quantity $B$ defined as:
\begin{equation}
B = \left(1+{2\alpha\over\beta}t^2 {1+At^2\over 1+At^2+A^2t^4 }\right),
\end{equation}
so that $H_0 = g^3(\beta^2/2\alpha)\mbox{log}B$.

The derivatives with respect to $r_s$:
\begin{equation}
{\partial H_0(r_s,t,\zeta)\over\partial r_s}=
   g^3{\beta^2\over 2\alpha}
      {1 \over B}
       {\partial B \over\partial A} 
           {\partial A\over\partial r_s}~,\qquad
\end{equation}
where
\begin{equation}
{\partial B\over\partial A} =
   -{2\alpha\over\beta}A t^6 
           {2+At^2\over (1+At^2+A^2t^4)^2}~\qquad
{\partial A\over\partial r_s}=
   -{4\alpha^2\over g^3\beta^3} 
    {e^{-2\alpha\epsilon_c(r_s,\zeta)/(g^3\beta^2)}
     \over 
     \left(e^{-2\alpha\epsilon_c(r_s,\zeta)/(g^3\beta^2)}-1\right)^2}
    {\partial \epsilon_c(r_s,\zeta)\over\partial r_s}
\end{equation}
\begin{equation}
{\partial H_1(r_s,t,\zeta)\over\partial r_s}=
 \nu {d C_c(r_s)\over d r_s}g^3 t^2 
      e^{-a_4 (k^2_s/k^2_F)g^4 t^2} - d_0 a_4 g^4 t^2 H_1
\end{equation}
where $d_0=0.663436444$ and
\begin{equation}
{d C_c(r_s)\over d r_s}=
 {(c_2 +2c_3 r_s)(1+c_4 r_s+c_5 r_s^2+c_6 r_s^3)-
  (c_1+c_2r_s +c_3 r_s^2)(c_4 +2c_5 r_s+3c_6 r_s^2)
  \over \left(1+c_4 r_s+c_5 r_s^2+c_6 r_s^3\right)^2}.
\end{equation}

The derivatives with respect to $\zeta$:
\begin{equation}
{\partial H_0(r_s,t,\zeta)\over\partial\zeta}=
 {3\over g} H_0{\partial g\over\partial\zeta} +
  g^3{\beta^2\over 2\alpha}
      {1 \over B}
       {\partial B\over\partial A} 
           {\partial A\over\partial\zeta}.
\end{equation}
where
\begin{equation}
{\partial A\over\partial\zeta}=
    -{4\alpha^2\over g^3\beta^3} 
    {e^{-2\alpha\epsilon_c(r_s,\zeta)/(g^3\beta^2)}
     \over 
     \left(e^{-2\alpha\epsilon_c(r_s,\zeta)/(g^3\beta^2)}-1\right)^2}
   \left({\partial \epsilon_c(r_s,\zeta)\over\partial\zeta}
         -{3\epsilon_c(r_s,\zeta)\over g}{\partial g\over\partial\zeta}
   \right)
\end{equation}
\begin{equation}
{\partial H_1(r_s,t,\zeta)\over\partial r_s}=
     \left( {3\over g}- 4a_4 (k^2_s/k^2_F) g^3 t^2\right)
     H_1 {\partial g\over\partial\zeta}.
\end{equation}

The derivatives with respect to $t$:
\begin{equation}
{\partial H_0(r_s,t,\zeta)\over\partial\zeta}=
   g^3{\beta^2\over 2\alpha}
      {1 \over B}
       {\partial B\over\partial t} 
\end{equation}
where
\begin{equation}
{\partial B\over\partial t} =
   {4\alpha\over\beta}t 
           {1+At^2\over 1+At^2+A^2t^4}+
   {2\alpha\over\beta}t^2
           {(1+2At)(1+At^2+A^2t^4)-(1+At^2)(2At+4A^2t^3)
            \over \left(1+At^2+A^2t^4\right)^2}
\end{equation}
\begin{equation}
{\partial H_1(r_s,t,\zeta)\over\partial r_s}=
     \left( {2\over t}- 2a_4 (k^2_s/k^2_F) g^4 t\right) H_1 .
\end{equation}

\end{document}