Introduction
In the PAW method, the all-electron (AE) wavefunction (AEWFC) $\psi_{n\mathbf{k}}$ is related to the pseudo-wavefucntion (PSWFC) $\tilde\psi_{n\mathbf{k}}$ by means of a linear transformation 1
\[\begin{equation} |\psi_{n\mathbf{k}}\rangle = |\tilde\psi_{n\mathbf{k}} \rangle + \sum_j [ |\phi_j\rangle - |\tilde\phi_j\rangle ] \cdot \langle \tilde{p}_j | \tilde\psi_{n\mathbf{k}} \rangle \label{eq:paw_wfc} \end{equation}\]Note that both AEWFC and PSWFC are defined on the Born-von Kármán (BvK) supercell rather than the unit cell, therefore, the index j in Eq.\eqref{eq:paw_wfc} depends not only on the atom index in the unit cell but also on the index of the unit cell, i.e.
\[\begin{equation} \sum_j \quad\Longrightarrow\quad \sum_{\mathbf{R}_L} \sum_{\boldsymbol\tau_\alpha} \sum_k \end{equation}\]The forms of the projector function $\tilde{p}$, and the partial waves $\phi$ and $\tilde\phi$, however, are only determined by the kind of the atom except that the centers depend on the cell and atom indices.
AEWFC and PSWFC
The AEWFC $\psi_{n\mathbf{k}}$ and the PSWFC $\tilde\psi_{n\mathbf{k}}$ both take the form of a Bloch wave, i.e.
\[\begin{align} |\psi_{n\mathbf{k}}\rangle &= e^{i\mathbf{k}\cdot\mathbf{r}} \cdot |u_{n\mathbf{k}}\rangle \\[6pt] |\tilde\psi_{n\mathbf{k}}\rangle &= e^{i\mathbf{k}\cdot\mathbf{r}} \cdot |\tilde{u}_{n\mathbf{k}}\rangle \end{align}\]where $u_{n\mathbf{k}}$ and $\tilde{u}_{n\mathbf{k}}$ are cell-periodic functions.
In VASP
, cell-periodic part of the PSWFCs, i.e.
$|\tilde{u}_{n\mathbf{k}}\rangle$, are the ones included in the WAVECAR
file,
only that the plane-waves coefficients $C_{n\mathbf{k}}(\mathbf{G})$ are
written instead.
where $N_\mathrm{fft}$ is the number of FFT grid points, $V$ is the cell volume, and the plane-waves are limited within the sphere defined by
\[\begin{equation} \frac{\hbar^2|\mathbf{G} + \mathbf{k}|^2}{2m_e} < E_c \label{eq:pw_encut} \end{equation}\]The relation between the norm of $\tilde{u}_{n\mathbf{k}}(\mathbf{r})$ and $C_{n\mathbf{k}}(\mathbf{G})$ are given by the following equation
\[\begin{equation} \sum_{\mathbf{G}} |C_{n\mathbf{k}}(\mathbf{G})|^2 = \frac{V}{N_\mathrm{fft}} \sum_{\mathbf{r}} |\tilde{u}_{n\mathbf{k}}(\mathbf{r})|^2 \approx 1 \end{equation}\]Note that the norm of the PSWFC can be larger or smaller than one, hence we use “$\approx$” instead of equal sign in the above equation.
To obtain the real-space representation of $\tilde{u}_{n\mathbf{k}}$, one need to perform Fourier transform (FT) on the plane-wave coefficients, which can be easily realized utilizing the VaspBandUnfolding package.
The real-space representation of $u_{n\mathbf{k}}$ is not as readily obtainable, since it contains the on-site contributions that are defined on multiple radial logarithmic grids with different centers, which are not compatible with the uniform grid used to represent $\tilde{u}_{n\mathbf{k}}$.
Projector Function
The projector functions $\tilde{p}_j$ depend on the distance from the center of the PAW sphere on which they are localized,
\[\begin{equation} \langle \mathbf{r} | \tilde{p}_j \rangle = \tilde{p}_j(\mathbf{r} - \mathbf{R}_L - \boldsymbol{\tau}_j); \qquad | \mathbf{r} - \mathbf{R}_L - \boldsymbol{\tau}_j | < r_j^c \end{equation}\]where $r_j^c$ is the PAW sphere cutoff radius, $\mathbf{R}_L$ is the cell vector of the L-th unit cell and $\boldsymbol\tau_j$ is the atom position within the cell.
Note in passing that $r_j^c$ in
VASP
does not necessarily satisfy the non-overlaping condition. For examples, $r_O^c = 0.82 \,\mathring{A}$ and $r_H^c = 0.59 \,\mathring{A}$ while a typical O-H bond length is about $0.97 \,\mathring{A}$.
For simplicity, we only consider the projector functions located within the first unit cell, i.e. $\mathbf{R}_0 = (0,0,0)$. The form of the projector function is composed of a radial part and an angular part, i.e.
\[\begin{equation} \tilde{p}_j(\mathbf{r} - \boldsymbol\tau) = f_j(|\mathbf{r} - \boldsymbol\tau|) \cdot Y_{lm}(\widehat{\mathbf{r} - \boldsymbol\tau}) \end{equation}\]where REAL FUNCTION $f_j(r)$ and real spherical harmonics $Y_{lm}$ are used in order to save computational cost. The Fourier transform of the projector function, as is shown in a previous post,2 can also be written as a radial function multiplied by a real spherical harmonics with the same $(l,m)$ quantum numbers, i.e.
\[\begin{align} {\cal F}[\tilde{p}_j(\mathbf{r} - \boldsymbol\tau)](\mathbf{k}) &= \frac{1}{\sqrt{V}} \int_V\, \tilde{p}_j(\mathbf{r} - \boldsymbol\tau) \, e^{-i\mathbf{k}\cdot\mathbf{r}} \, \mathrm{d}\mathbf{r} \\[6pt] &= \frac{i^{-\ell}}{\sqrt{V}}\cdot g_j(k) \cdot Y_{lm}(\hat{\mathbf{k}}) \cdot e^{-i\mathbf{k}\cdot\boldsymbol\tau} \end{align}\]where the phase term is due to the shifting property of FT. $f_j(r)$ and $g_j(k)$ are linked by spherical Bessel transform,3
\[\begin{equation} g_j(k) = 4\pi \int_0^\infty f_j(r) \,j_l(kr) r^2 \,\mathrm{d} r \end{equation}\]Obviously, $g_j(k)$ is also a REAL FUNCTION by definition.
In VASP
, the information of the projector functions are included in the
POTCAR
file, however, only radial parts of the real- and momentum-space
projector functions, i.e. $f_j(r)$ and $g_j(k)$, are are stored. As a result,
one must multiply the radial function with the corresponding real spherical
harmonics $Y_{lm}$ in order to get the complete projector function.
To evaluate the inner-product of the projector function and PSWFC, an extra
k-dependent projector function $\tilde{p}_{j\mathbf{k}}$ is introduced in
VASP
4
As a result, the inner product between the projector and PSWFC can be written as
\[\begin{align} \beta_{n\mathbf{k}}^j &= \langle \tilde{p}_j | \tilde\psi_{n\mathbf{k}} \rangle \\[6pt] &= e^{i\mathbf{k}\cdot\boldsymbol{\tau}_j} \cdot \langle \tilde{p}_{j\mathbf{k}} | \tilde{u}_{n\mathbf{k}} \rangle \\[6pt] &= e^{i\mathbf{k}\cdot\boldsymbol{\tau}_j} \cdot \beta_{n\mathbf{k}}^{j\mathbf{k}} \label{eq:kdep_proj_inner} \end{align}\]Note that in practical implementation,
VASP
DOES NOT include the phase term $e^{i\mathbf{k}\cdot\boldsymbol{\tau}_j}$ in the actual calculation, since it does not change the occupation matrix or energy, as is explained in thenonl.F
ofVASP
source file.
One can clearly see from Eq.\eqref{eq:kdep_proj_inner} that for the same projector function located at another unit cell $\mathbf{R}_L$, the inner product is simply right-hand-side of Eq.\eqref{eq:kdep_proj_inner} multiplied by $e^{i\mathbf{k}\cdot\mathbf{R}_L}$.
With this definition, Eq.\eqref{eq:kdep_proj_inner} can be evaluated in two different approaches:
- momentum-space method
In the first step, one has to find out all the plane-waves within the sphere defined by Eq.\eqref{eq:pw_encut}. Secondly, $g_j$ values on the relevant $|\mathbf{G} + \mathbf{k}|$ are obtained by Spline-interpolation and multiplied with the corresponding real spherical harmonics $Y_{lm}$. The sum over $\mathbf{G}$ can then readily obtained with the knowledge of the plane-wave coefficients.
Note that this is what
VASP
does by settingLREAL = .F.
in theINCAR
.
- real-space method
First, the grid points within the PAW sphere $|\mathbf{r} - \boldsymbol\tau_j| < r_j^c$ are identified. Secondly, the values of $f_j$ on these grid points are obtained by Spline-interpolation and multiplied with the corresponding real spherical harmonics $Y_{lm}$. Finally, muliplying the phase term and $\tilde{u}_{n\mathbf{k}}(\mathbf{r})$, the integral becomes a sum over the selected grid points.
Note that this is what
VASP
does by settingLREAL = .T.
in theINCAR
.
Partial Waves
The AE/PS partial waves $\phi_j$ and $\tilde\phi_j$ are defined on the radial logarithmic grid with center at $\boldsymbol\tau_j$ (in the first unit cell), i.e.
\[\begin{align} \langle \mathbf{r} | \phi_j \rangle &= \phi_j(\mathbf{r} - \boldsymbol{\tau}_j) \\[6pt] \langle \mathbf{r} | \tilde{\phi}_j \rangle &= \tilde{\phi}_j(\mathbf{r} - \boldsymbol{\tau}_j) \end{align}\]The radial logarithmic grid is defined as
\[r_{n} = r_0 \cdot e^{n\cdot H};\qquad n=0,1,2,\ldots,N-1\]or equivalently
\[r_{n+1} = r_n \cdot e^{H};\qquad n=0,1,2,\ldots,N-1\]where $r_0$ is the minimum value of $|\mathbf{r} - \boldsymbol\tau_j|$ and is by definition not zero, $N$ is the number of radial grid points.
Note that $r_{N-1}$ is not necessarily equal to the PAW sphere cutoff radius $r_j^c$, as can be seen in Figure 2. However, $\phi_j$ and $\tilde\phi_j$ are indential beyond $|\mathbf{r} - \boldsymbol\tau_j| = r_j^c$, which makes the difference of the two zero in Eq.\eqref{eq:paw_wfc} beyond $r_j^c$.
The AE/PS partial waves are again composed of radial and angular parts,
\[\begin{align} \phi_j(\mathbf{r} - \boldsymbol{\tau}_j) = \frac{ R_j(|\mathbf{r} - \boldsymbol\tau_j|) }{ |\mathbf{r} - \boldsymbol\tau_j| } \cdot Y_{lm}(\widehat{\mathbf{r} - \boldsymbol\tau_j}) \\[6pt] \tilde\phi_j(\mathbf{r} - \boldsymbol{\tau}_j) = \frac{ \tilde{R}_j(|\mathbf{r} - \boldsymbol\tau_j|) }{ |\mathbf{r} - \boldsymbol\tau_j| } \cdot Y_{lm}(\widehat{\mathbf{r} - \boldsymbol\tau_j}) \end{align}\]IT MUST BE STRESSED that $R_j(r)$ and $\tilde{R}_j(r)$ are the ones stored in the
POTCAR
file and are also the ones shown in Figure 2.
A few notes on the construction of AE/PS partial waves
-
The AE partial waves are chosen to be solution to Schrödinger equation for a fixed quantum number l and chosen reference energy. 5
-
The PS partial waves are, by definition, indential to the AE partial waves outside th augmentation sphere. Inside the sphere, the PS partial waves are chosen to be a smooth continuation of the AE ones.
\[\begin{equation} \frac{\tilde{R}_j(r)}{r} = \begin{cases} \displaystyle \sum_{i=1}^2 \alpha_i j_l(q_k r), & r \lt r_c \\[6pt] \displaystyle \frac{R_j(r)}{r}, & r \gt r_c \end{cases} \end{equation}\]VASP
expands the PS partial waves within the sphere in terms of a linear combination of spherical Bessel functions 5 6where the coefficients $\alpha_i$ and $q_i$ are determined so that the PS partial waves are two times continuously differentiable at $r_c$.
Alternatively, one can use polynomials of r with even powers 7 (Troullier-Martins pseudopotential 8) or Gaussian in the expansion.
Constructing AEWFC
As has been stated, the AEWFC $\psi_{n\mathbf{k}}$ is a Bloch wave and is defined on the BvK supercell. As a result, one can only construct the AEWFC within a few unit cells if one wishes to get the real-space representation of the AEWFC.
Alternatively, one can construct the cell-periodic part of the AEWFC $u_{n\mathbf{k}}(\mathbf{r})$
\[\begin{align} |u_{n\mathbf{k}}\rangle &= |\tilde{u}_{n\mathbf{k}} \rangle + e^{-i\mathbf{k}\cdot\mathbf{r}} \sum_j \left[ |\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle - |\tilde\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle \right] \cdot \langle \tilde{p}_j(\mathbf{r} - \boldsymbol\tau_j) | \tilde\psi_{n\mathbf{k}} \rangle \nonumber \\[9pt] &= |\tilde{u}_{n\mathbf{k}} \rangle + \sum_j e^{-i\mathbf{k}\cdot(\mathbf{r} - \boldsymbol\tau_j)} \left[ |\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle - |\tilde\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle \right] \cdot \beta_{n\mathbf{k}}^{j\mathbf{k}} \label{eq:paw_periodic_wfc} \end{align}\]Similar to the evaluation of Eq.\eqref{eq:kdep_proj_inner}, Eq.\eqref{eq:paw_periodic_wfc} can also be performed in real- and momentum-space.
IT MUST BE STRESSED that although the all-electron quantities can be obtained in principle, they are NEVER handled directly in practical calculations. The essence of the PAW method lies in the fact that there are NO CROSS TERMS between quantities on the regular (coarse) uniform grid and the quantities on the radial (fine) grid, and as a result the two kinds of quantities can be treated separately.
First of all, set the size of the new uniform grid that is used to represent $u_{n\mathbf{k}}$. This can be done by introducting a new energy cutoff $E_{\mathrm{ae}}$, which must be larger than the energy cutoff $E_c$ that determines the size of the uniform grid used to represent $\tilde{u}_{n\mathbf{k}}$. Obviously, the new and old grid size are related by
\[\begin{equation} \begin{pmatrix} n_x^{\mathrm{ae}},& n_y^{\mathrm{ae}},& n_z^{\mathrm{ae}} \end{pmatrix} = \left\lceil \sqrt{\frac{E_{\mathrm{ae}}}{E_{\mathrm{c}}}} \cdot \begin{pmatrix} n_x^{\mathrm{c}},& n_y^{\mathrm{c}},& n_z^{\mathrm{c}} \end{pmatrix} \right\rceil \end{equation}\]where $\lceil\ldots\rceil$ denotes a integer ceiling operation.
Secondly, Fourier transform the plane-wave ceofficients $C_{n\mathbf{k}}(\mathbf{G})$ to get $\tilde{u}_{n\mathbf{k}}$ on the new grid.
The rest is to get the real-space representation of the on-site contribution (the summation term in Eq.\eqref{eq:paw_periodic_wfc}), which can be evaluated by two approaches.
real-space approach
The on-site term in real space is simply
\[\begin{align} &\phantom{={}}\sum_j \beta_{n\mathbf{k}}^{j\mathbf{k}} \cdot e^{-i\mathbf{k}\cdot(\mathbf{r} - \boldsymbol\tau_j)} \cdot \left[ \phi_j(\mathbf{r} - \boldsymbol\tau_j) - \tilde\phi_j(\mathbf{r} - \boldsymbol\tau_j) \right] \nonumber \\[6pt] &=\sum_j \beta_{n\mathbf{k}}^{j\mathbf{k}} \cdot e^{-i\mathbf{k}\cdot(\mathbf{r} - \boldsymbol\tau_j)} \cdot \frac{ \left[ R_j(\mathbf{r} - \boldsymbol\tau_j) - \tilde{R}_j(\mathbf{r} - \boldsymbol\tau_j) \right] }{ |\mathbf{r} - \boldsymbol\tau_j| } \cdot Y_{lm}(\widehat{\mathbf{r} - \boldsymbol\tau_j}) \end{align}\]This can be done similar to the the real-space evaluation of $\beta_{n\mathbf{k}}^{j\mathbf{k}}$.
momentum-space approach
Let’s now inspect the plane-wave expansion of the on-site terms
\[\begin{align} S_j(\mathbf{G}) &= \beta_{n\mathbf{k}}^{j\mathbf{k}} \cdot \langle \mathbf{G} | e^{-i\mathbf{k}\cdot(\mathbf{r} - \boldsymbol\tau_j)} \left[ |\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle - |\tilde\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle \right] \nonumber \\[9pt] &= \beta_{n\mathbf{k}}^{j\mathbf{k}} \, e^{i\mathbf{k}\cdot\boldsymbol\tau_j} \cdot \sum_{\mathbf{G}'} \langle \mathbf{G} | e^{-i\mathbf{k}\cdot\mathbf{r}} | \mathbf{G}' \rangle \langle \mathbf{G}' | \left[ |\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle - |\tilde\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle \right] \nonumber \\[6pt] &= \beta_{n\mathbf{k}}^{j\mathbf{k}} \, e^{i\mathbf{k}\cdot\boldsymbol\tau_j} \cdot \sum_{\mathbf{G}'} \delta_{\mathbf{G} + \mathbf{k}, \mathbf{G}'} \cdot \langle \mathbf{G}' | \left[ |\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle - |\tilde\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle \right] \nonumber \\[6pt] &= \beta_{n\mathbf{k}}^{j\mathbf{k}} \, e^{i\mathbf{k}\cdot\boldsymbol\tau_j} \cdot \langle \mathbf{G} + \mathbf{k} | \left[ |\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle - |\tilde\phi_j(\mathbf{r} - \boldsymbol\tau_j)\rangle \right] \nonumber \\[6pt] &= \frac{i^{-l}}{\sqrt{V}}\, \beta_{n\mathbf{k}}^{j\mathbf{k}} \, e^{i\mathbf{k}\cdot\boldsymbol\tau_j} \cdot \left[ {\cal G}_j(|\mathbf{G} + \mathbf{k}|) - \tilde{\cal G}_j(|\mathbf{G} + \mathbf{k}|) \right] \cdot Y_{lm}(\widehat{\mathbf{G} + \mathbf{k}}) \cdot e^{-i(\mathbf{G}+\mathbf{k})\cdot\boldsymbol\tau_j} \nonumber \\[6pt] &= \frac{i^{-l}}{\sqrt{V}}\, \beta_{n\mathbf{k}}^{j\mathbf{k}} \cdot \left[ {\cal G}_j(|\mathbf{G} + \mathbf{k}|) - \tilde{\cal G}_j(|\mathbf{G} + \mathbf{k}|) \right] \cdot Y_{lm}(\widehat{\mathbf{G} + \mathbf{k}}) \cdot e^{-i\mathbf{G}\cdot\boldsymbol\tau_j} \end{align}\]where $\cal G$ is the spherical Bessel transform of the partical waves
\[\begin{align} {\cal G}_j(k) &= 4\pi \int_0^\infty \frac{R_j(r)}{r} \,j_l(kr) r^2 \,\mathrm{d} r \\[6pt] \tilde{\cal G}_j(k) &= 4\pi \int_0^\infty \frac{\tilde{R}_j(r)}{r} \,j_l(kr) r^2 \,\mathrm{d} r \end{align}\]Note that the plane-waves are now determined by
\[\begin{equation} \frac{\hbar^2|\mathbf{G} + \mathbf{k}|^2}{2m_e} < E_{\mathrm{ae}} \label{eq:pw_encut_ae} \end{equation}\]The real-space on-site term can be obtained by Fourier transform of $\sum_j S_j(\mathbf{G})$.
AEWFC & PSWFC of CO2 Orbital
Here, we choose CO2 as an example to investigate the AEWFC and PSWFC.
In principle, the norm of AEWFC is unity and that of the PSWFC is different from one according to Eq.\eqref{eq:paw_wfc}. Therefore, to best show the difference between the two wavefunctions, we first interate through the bands and find out the orbital with the maximum contribution from the PAW core region.
1
2
3
4
5
6
7
8
9
10
11
12
13
import numpy as np
from vaspwfc import vaspwfc
# the pseudo-wavefunction
ps_wfc = vaspwfc('WAVECAR', lgamma=True)
# Find out the band that has the most core contributions
for iband in range(ps_wfc._nbands):
# plane-wave coefficients of PSWFC
cg = ps_wfc.readBandCoeff(iband=iband+1)
# norm of the PSWFC
ps_norm = np.sum(cg.conj() * cg).real
print(f"#band: {iband+1:3d} -> {1 - ps_norm: 8.4f}")
The output is shown below, where one can see that band 7 and 8 possess the maximum contributions from the core region. One can also see that the contribution can be both positive and negative, which means that the norm of the PSWFC can be smaller (positive output) or larger (negative output) than unity.
1
2
3
4
5
6
7
8
9
10
11
12
#band: 1 -> -0.1033
#band: 2 -> -0.0744
#band: 3 -> 0.0302
#band: 4 -> 0.1362
#band: 5 -> 0.1362
#band: 6 -> 0.1391
#band: 7 -> 0.1993
#band: 8 -> 0.1993
#band: 9 -> -0.0010
#band: 10 -> 0.1494
#band: 11 -> 0.1494
#band: 12 -> 0.0067
We will choose band 8, which is also the HOMO orbital, to proceed. Below,
aecut=-25
means that $E_\mathrm{ae} = 25\times E_\mathrm{c}$.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
from vaspwfc import vaspwfc
from aewfc import vasp_ae_wfc
# the pseudo-wavefunction
ps_wfc = vaspwfc('WAVECAR', lgamma=True)
# the all-electron wavefunction
ae_wfc = vasp_ae_wfc(ps_wfc, aecut=-25)
which_band = 8
phi_ae, phi_core_ae, phi_core_ps = ae_wfc.get_ae_wfc(
iband=which_band, lcore=True
)
phi_ps = ps_wfc.get_ps_wfc(
iband=which_band,
norm=False, ngrid=ae_wfc._aegrid
)
Below, we show the difference between PSWFC and AEWFC of CO2 HOMO.
-
PSWFC, AEWFC and Core AE/PS Partial WFC
-
Wavefunctions $\psi(r,z)$ along the bond direction $z$ at different vertical distances $r$.
Reference
-
pySBT: A python implementation of spherical Bessel transform (SBT) in O(Nlog(N)) time. ↩
-
Linear optical properties in the projector-augmented wave methodology ↩
-
From ultrasoft pseudopotentials to the projector augmented-wave method ↩ ↩2
-
In another VASP lecture note, 3-4 spherical Bessel functions are used in the expansion. ↩
-
As odd powers in r results in a kink in the functions at the origin, i.e. that the first derivatives are not defined at this point. ↩