PAW All-Electron Wavefunction in VASP

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 ¹

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.

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.

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

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

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,

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.

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,² can also be written as a radial function multiplied by a real spherical harmonics with the same $(l,m)$ quantum numbers, i.e.

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,³

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⁴

As a result, the inner product between the projector and PSWFC can be written as

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 the nonl.F of VASP 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 setting LREAL = .F. in the INCAR.

• 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 setting LREAL = .T. in the INCAR.

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.

The radial logarithmic grid is defined as

or equivalently

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,

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. ⁵

• 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. VASP expands the PS partial waves within the sphere in terms of a linear combination of spherical Bessel functions ^{5 6}

  \[\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}\]

  where 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 ⁷ (Troullier-Martins pseudopotential ⁸) 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})$

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

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

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

where $\cal G$ is the spherical Bessel transform of the partical waves

Note that the plane-waves are now determined by

The real-space on-site term can be obtained by Fourier transform of $\sum_j S_j(\mathbf{G})$.

AEWFC & PSWFC of CO₂ Orbital

Here, we choose CO₂ 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 CO₂ HOMO.

• PSWFC, AEWFC and Core AE/PS Partial WFC

  

• Wavefunctions $\psi(r,z)$ along the bond direction $z$ at different vertical distances $r$.

  

Reference