# Developers guide to GPAW¶

XXX Update page to new GPAW style (after guc merge) and mention NewLFCs.

This page goes through the most important equations of a PAW calculation and has references to the code. It is a good idea to have the big picture in front of you when reading this page.

• Initial wave functions and densities (todo)

• Finding the ground state (todo)

## Wave functions¶

The central quantities in a PAW calculation are the pseudo wave-functions, $$\tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r})$$, from which the all-electron wave functions can be obtained:

$\psi_{\sigma\mathbf{k}n}(\mathbf{r}) = \tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r}) + \sum_a \sum_i [\phi_i^a(\mathbf{r} - \mathbf{R}^a) - \tilde{\phi}_i^a(\mathbf{r} - \mathbf{R}^a)] \langle\tilde{p}_i^a | \tilde{\psi}_{\sigma\mathbf{k}n} \rangle,$

where

$\langle\tilde{p}_i^a | \tilde{\psi}_{\sigma\mathbf{k}n} \rangle = \int d\mathbf{r} \tilde{p}_i^a(\mathbf{r} - \mathbf{R}^a) \tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r}).$

Here, $$a$$ is the atom number, $$\mathbf{R}^a$$ is the position of atom number $$a$$ and $$\tilde{p}_i^a$$, $$\tilde{\phi}_i^a$$ and $$\phi_i^a$$ are the projector functions, pseudo partial waves, and all-electron partial waves respectively, of the atoms.

See Naming convention for arrays for more information on the naming of arrays. Note that spos_c gives the position of the atom in scaled coordinates in the range [0:1[ (relative to the unit cell).

Note, that in the code, i refers to $$n$$, $$\ell$$ and $$m$$ quantum numbers, and j refers to $$n$$ and $$\ell$$ only (see Naming convention for arrays). So, to put an atom-centered function like $$\tilde{p}_{n\ell m}^a(\mathbf{r})$$ on the 3D grid, you need both the radial part $$\tilde{p}_{n\ell}^a(r)$$ (one of the splines in paw.wfs.setups[a].pt_j) and a spherical harmonics $$Y_{\ell m}(\theta,\phi)$$. Putting radial functions times spherical harmonics on a grid is done by the gpaw.lfc.LocalizedFunctionsCollection class. See also gpaw.setup.Setup and gpaw.spline.Spline.

The wave-functions are othonormalized such that the pseudo wave-functions obey the following orthogonality requirements:

$\langle \psi_{\sigma\mathbf{k}n} | \psi_{\sigma\mathbf{k}m} \rangle = \langle \tilde{\psi}_{\sigma\mathbf{k}n} | \hat{O} | \tilde{\psi}_{\sigma\mathbf{k}m} \rangle = \delta_{nm},$

, where $$\hat{O}$$ is the overlap operator in the PAW formalism. Refer to Orthogonalizing the wave functions for details.

## Overlaps¶

The overlap operator is defined in terms of the PAW overlap corrections:

$\hat{O} = 1 + \sum_a \sum_{i_1 i_2} |\tilde{p}_{i_1}^a\rangle \Delta O_{i_1 i_2}^a \langle\tilde{p}_{i_2}^a|.$

The constants $$\Delta O_{i_1 i_2}^a$$ are found in paw.wfs.setups[a].dO_ii (ndarray). XXX Someone should rename dO_ii to dS_ii or $$\hat{S}$$ to $$\hat{O}$$.

$\Delta O_{i_1 i_2}^a = \int d\mathbf{r} [\phi_{i_1}^a(\mathbf{r})\phi_{i_2}^a(\mathbf{r}) - \tilde{\phi}_{i_1}^a(\mathbf{r})\tilde{\phi}_{i_2}^a(\mathbf{r})].$

An approximate inverse overlap operator is similarly defined by:

$\hat{O}^{\;-1}_\mathrm{approx.} = 1 + \sum_a \sum_{i_1 i_2} |\tilde{p}_{i_1}^a\rangle \Delta C_{i_1 i_2}^a \langle\tilde{p}_{i_2}^a|.$

The inverse overlap coefficients $$\Delta C_{i_1 i_2}^a$$ are found in setup.dC_ii (ndarray) and are solutions to the system of linear equations:

$\Delta C_{i_1 i_2}^a + \Delta O_{i_1 i_2}^a + \sum_{i_3 i_4} \Delta C_{i_1 i_3}^a B_{i_3 i_4}^a \Delta O_{i_4 i_2}^a = 0 \qquad ,\forall i_1,i_2$

, such that $$\hat{O}^{\;-1}_\mathrm{approx.}\hat{O} = \hat{I}$$ provided $$\langle\tilde{p}_{i_1}^a|\tilde{p}_{i_2}^{a'}\rangle = \delta_{a a'} \langle\tilde{p}_{i_1}^a|\tilde{p}_{i_2}^{a}\rangle$$. These projector overlaps $$B_{i_1 i_2}^a = \langle\tilde{p}_{i_1}^a|\tilde{p}_{i_2}^{a}\rangle$$ are likewise found in setup.B_ii.

## Densities¶

From the pseudo wave-functions, the pseudo electron spin-densities can be constructed (see here):

$\tilde{n}_\sigma(\mathbf{r}) = \frac{1}{N_s} \sum_{s=1}^{N_s} \hat{S}_s \left [ \sum_{n\mathbf{k}} f_{n\mathbf{k}\sigma} |\tilde{\psi}_{n\mathbf{k}\sigma}(\mathbf{r})|^2 + \frac{1}{2} \sum_a \tilde{n}_c^a(|\mathbf{r}-\mathbf{R}^a|) \right ].$

Here, $$\hat{S}_s$$ is one of the $$N_s$$ symmetry operators of the system (see gpaw.symmetry.Symmetry), $$f_{n\mathbf{k}\sigma}$$ are the occupation numbers (adding up to the number of valence elctrons), and $$\tilde{n}_c^a(r)$$ is the pseudo core density for atom number $$a$$.

The all-electron spin-densities are given as:

$n_\sigma(\mathbf{r}) = \tilde{n}_\sigma(\mathbf{r}) + \sum_a [n_\sigma^a(\mathbf{r} - \mathbf{R}^a) - \tilde{n}_\sigma^a(\mathbf{r} - \mathbf{R}^a)],$

where

$n_\sigma^a(\mathbf{r}) = \sum_{i_1 i_2} D_{\sigma i_1 i_2}^a \phi_{i_1}^a(\mathbf{r})\phi_{i_2}^a(\mathbf{r}) + \frac{1}{2} n_c^a(r),$
$\tilde{n}_\sigma^a(\mathbf{r}) = \sum_{i_1 i_2} D_{\sigma i_1 i_2}^a \tilde{\phi}_{i_1}^a(\mathbf{r})\tilde{\phi}_{i_2}^a(\mathbf{r}) + \frac{1}{2} \tilde{n}_c^a(r),$

are atom centered expansions, and

$D_{\sigma i_1 i_2}^a = \sum_{n\mathbf{k}} \langle \tilde{\psi}_{\sigma\mathbf{k}n} | \tilde{p}_{i_1}^a \rangle f_{n\mathbf{k}\sigma} \langle \tilde{p}_{i_2}^a | \tilde{\psi}_{\sigma\mathbf{k}n} \rangle$

is an atomic spin-density matrix, which must be symmetrized the same way as the pseudo electron spin-densities.

 formula object type $$\hat{S}_s$$ paw.wfs.symmetry gpaw.symmetry.Symmetry $$\tilde{n}_\sigma$$ paw.density.nt_sG and paw.density.nt_sg ndarray $$\tilde{n}=\sum_\sigma\tilde{n}_\sigma$$ paw.density.nt_g ndarray $$\tilde{n}_c^a(r)$$ paw.wfs.setups[a].nct gpaw.spline.Spline $$\tilde{n}_c^a(\mathbf{r}-\mathbf{R}^a)$$ paw.density.nct gpaw.lfc.LocalizedFunctionsCollection $$f_{\sigma\mathbf{k}n}$$ paw.wfs.kpt_u[u].f_n ndarray $$D_{\sigma i_1 i_2}^a$$ paw.density.D_asp[a] ndarray

From the all-electron and pseudo electron densities we can now construct corresponding total all-electron and pseudo charge densities:

$\rho(\mathbf{r}) = \sum_\sigma n_\sigma(\mathbf{r}) + \sum_a Z^a(\mathbf{r} - \mathbf{R}^a),$
$\tilde{\rho}(\mathbf{r}) = \sum_\sigma \tilde{n}_\sigma(\mathbf{r}) + \sum_a \tilde{Z}^a(\mathbf{r} - \mathbf{R}^a).$

If $$\mathbb{Z}^a$$ is the atomic number of atom number $$a$$, then $$Z^a(\mathbf{r})=-\mathbb{Z}^a\delta(\mathbf{r})$$ (we count the electrons as positive charge and the protons as negative charge). The compensation charges are given as:

$\tilde{Z}^a(\mathbf{r}) = \sum_{\ell=0}^{\ell_{\text{max}}} \sum_{m=-\ell}^\ell Q_{\ell m}^a \hat{g}_{\ell m}^a(\mathbf{r}) = \sum_{\ell=0}^{\ell_{\text{max}}} \sum_{m=-\ell}^\ell Q_{\ell m}^a \hat{g}_\ell^a(r) Y_{\ell m}(\theta,\phi),$

where $$\hat{g}_\ell^a(r)\propto r^\ell\exp(-\alpha^a r^2)$$ are Gaussians. The compensation charges should make sure that the two atom centered densities $$\rho^a=\sum_\sigma n_\sigma^a + Z^a$$ and $$\tilde{\rho}^a=\sum_\sigma \tilde{n}_\sigma^a + \tilde{Z}^a$$ have identical multipole expansions outside the augmentation sphere. This gives the following equation for $$Q_L^a$$:

$Q_L^a = \sum_{i_1 i_2} \Delta_{i_1 i_2 L}^a \sum_\sigma D_{\sigma i_1 i_2}^a + \Delta_0^a \delta_{\ell,0},$

where

$\Delta_{i_1 i_2 L}^a = \int d\mathbf{r} Y_L(\hat{\mathbf{r}}) r^\ell [\phi_{i_1}^a(\mathbf{r})\phi_{i_2}^a(\mathbf{r}) - \tilde{\phi}_{i_1}^a(\mathbf{r})\tilde{\phi}_{i_2}^a(\mathbf{r})],$
$\Delta_0^a = \int d\mathbf{r} Y_{00}(\hat{\mathbf{r}}) [-\mathbb{Z}^a \delta(\mathbf{r}) + n_c^a(\mathbf{r}) - \tilde{n}_c^a(\mathbf{r})].$
 formula object type $$\tilde{\rho}$$ paw.density.rhot_g ndarray $$\mathbb{Z}^a$$ setup.Z int $$\Delta_{i_1 i_2 L}^a$$ setup.Delta_pL ndarray $$\Delta_0^a$$ setup.Delta0 float $$\hat{g}_\ell^a(r)$$ setup.ghat_l List of gpaw.spline.Splines $$\hat{g}_L^a(\mathbf{r}-\mathbf{R}^a)$$ paw.density.ghat gpaw.lfc.LocalizedFunctionsCollection $$Q_L^a$$ paw.density.Q_aL[a] ndarray

## The total energy¶

The total PAW energy is composed of a smooth part evaluated using pseudo quantities on the 3D grid, plus corrections for each atom evaluated on radial grids inside the augmentation spheres: $$E=\tilde{E}+\sum_a(E^a - \tilde{E}^a)$$.

$\begin{split}\tilde{E} &= -\frac{1}{2} \sum_{\sigma\mathbf{k}n} f_{\sigma\mathbf{k}n} \int d\mathbf{r} \tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r}) \nabla^2 \tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r}) + \frac{1}{2}\int d\mathbf{r}d\mathbf{r}' \frac{\tilde{\rho}(\mathbf{r})\tilde{\rho}(\mathbf{r}')} {|\mathbf{r}-\mathbf{r}'|} \\ &\quad+ \sum_\sigma\sum_a\int d\mathbf{r}\tilde{n}_\sigma(\mathbf{r}) \bar{v}^a(|\mathbf{r}-\mathbf{R}^a|) + E_{\text{xc}}[\tilde{n}_\uparrow, \tilde{n}_\downarrow] % %.. math:: % \\ E^a &= -\frac{1}{2} 2\sum_i^{\text{core}} \int d\mathbf{r} \phi_i^a(\mathbf{r}) \nabla^2 \phi_i^a(\mathbf{r}) -\frac{1}{2} \sum_\sigma \sum_{i_1 i_2} D_{\sigma i_1 i_2}^a \int d\mathbf{r} \phi_{i_1}^a(\mathbf{r}) \nabla^2 \phi_{i_2}^a(\mathbf{r}) \\ &\quad+ \frac{1}{2}\int d\mathbf{r}d\mathbf{r}' \frac{\rho^a(\mathbf{r})\rho^a(\mathbf{r}')} {|\mathbf{r}-\mathbf{r}'|} + E_{\text{xc}}[n^a_\uparrow, n^a_\downarrow] % %.. math:: % \\ \tilde{E}^a &= -\frac{1}{2} \sum_\sigma\sum_{i_1 i_2} D_{\sigma i_1 i_2}^a \int d\mathbf{r} \tilde{\phi}_{i_1}^a(\mathbf{r}) \nabla^2 \tilde{\phi}_{i_2}^a(\mathbf{r}) + \frac{1}{2}\int d\mathbf{r}d\mathbf{r}' \frac{\tilde{\rho}^a(\mathbf{r})\tilde{\rho}^a(\mathbf{r}')} {|\mathbf{r}-\mathbf{r}'|} \\ &\quad+ \sum_\sigma \int d\mathbf{r}\tilde{n}^a_\sigma(\mathbf{r}) \bar{v}^a(r) + E_{\text{xc}}[\tilde{n}^a_\uparrow, \tilde{n}^a_\downarrow]\end{split}$

In the last two equations, the integrations are limited to inside the augmentation spheres only.

The electrostatic energy part of $$\tilde{E}$$ is calculated as $$\frac{1}{2}\int d\mathbf{r}\tilde{v}_H(\mathbf{r})\tilde{\rho}(\mathbf{r})$$, where the Hartree potential is found by solving Poissons equation: $$\nabla^2 \tilde{v}_H(\mathbf{r})=-4\pi\tilde{\rho}(\mathbf{r})$$ (see gpaw.poisson.FDPoissonSolver).