Posts PAW All-Electron Wavefunction in VASP
Post
Cancel

PAW All-Electron Wavefunction in VASP

Introduction

In the PAW method, the all-electron (AE) wavefunction (AEWFC) ψnk is related to the pseudo-wavefucntion (PSWFC) ψ~nk by means of a linear transformation 1

(1)|ψnk=|ψ~nk+j[|ϕj|ϕ~j]p~j|ψ~nk
Figure 1. Compositions of the all-electron wavefunction within the PAW method.

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.(1) depends not only on the atom index in the unit cell but also on the index of the unit cell, i.e.

(2)jRLταk

The forms of the projector function p~, and the partial waves ϕ and ϕ~, 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 ψnk and the PSWFC ψ~nk both take the form of a Bloch wave, i.e.

(3)|ψnk=eikr|unk(4)|ψ~nk=eikr|u~nk

where unk and u~nk are cell-periodic functions.

In VASP, cell-periodic part of the PSWFCs, i.e. |u~nk, are the ones included in the WAVECAR file, only that the plane-waves coefficients Cnk(G) are written instead.

Cnk(Gj)=u~nk(r)1VeiGrdr=ku~nk(rk)1VeiGjrkΔV(5)=VNfftku~nk(rk)1VeiGjrk(6)u~nk(rj)=kCnk(Gk)1VeiGkrj

where Nfft is the number of FFT grid points, V is the cell volume, and the plane-waves are limited within the sphere defined by

(7)2|G+k|22me<Ec

The relation between the norm of u~nk(r) and Cnk(G) are given by the following equation

(8)G|Cnk(G)|2=VNfftr|u~nk(r)|21

Note that the norm of the PSWFC can be larger or smaller than one, hence we use “” instead of equal sign in the above equation.

To obtain the real-space representation of u~nk, 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 unk 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 u~nk.

Projector Function

The projector functions p~j depend on the distance from the center of the PAW sphere on which they are localized,

(9)r|p~j=p~j(rRLτj);|rRLτj|<rjc

where rjc is the PAW sphere cutoff radius, RL is the cell vector of the L-th unit cell and τj is the atom position within the cell.

Note in passing that rjc in VASP does not necessarily satisfy the non-overlaping condition. For examples, rOc=0.82A˚ and rHc=0.59A˚ while a typical O-H bond length is about 0.97A˚.

For simplicity, we only consider the projector functions located within the first unit cell, i.e. R0=(0,0,0). The form of the projector function is composed of a radial part and an angular part, i.e.

(10)p~j(rτ)=fj(|rτ|)Ylm(rτ^)

where REAL FUNCTION fj(r) and real spherical harmonics Ylm 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.

(11)F[p~j(rτ)](k)=1VVp~j(rτ)eikrdr(12)=iVgj(k)Ylm(k^)eikτ

where the phase term is due to the shifting property of FT. fj(r) and gj(k) are linked by spherical Bessel transform,3

(13)gj(k)=4π0fj(r)jl(kr)r2dr

Obviously, gj(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. fj(r) and gj(k), are are stored. As a result, one must multiply the radial function with the corresponding real spherical harmonics Ylm in order to get the complete projector function.

Figure 2. Upper panels: radial parts of the projector functions in real (left) and momentum (right) space for Ti POTCAR. The vertical dashed lines in the left and right panels correspond to the real- and reciprocal-space cut, respectively. Lower panel: radial parts of the AE and PS partial waves on the logarithmic radial grid, shown with solid/dashed lines, respectively. The figure was generated with the potplot script.

To evaluate the inner-product of the projector function and PSWFC, an extra k-dependent projector function p~jk is introduced in VASP4

(14)|p~jk=eik(rτj)|p~j(rτj)

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

(15)βnkj=p~j|ψ~nk(16)=eikτjp~jk|u~nk(17)=eikτjβnkjk

Note that in practical implementation,VASP DOES NOT include the phase term eikτ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.(17) that for the same projector function located at another unit cell RL, the inner product is simply right-hand-side of Eq.(17) multiplied by eikRL.

With this definition, Eq.(17) can be evaluated in two different approaches:

  • momentum-space method
p~jk|u~nk=Gp~jk|GG|u~nk=Gp~j(rτj)|eik(rτj)|GG|u~nk=G,Gp~j(rτj)|GG|eik(rτj)|GG|u~nk=eikτjG,Gp~j(rτj)|GδG,G+kG|u~nk=eikτjGp~j(rτj)|G+kG|u~nk=eikτjVGilgj(|G+k|)Ylm(G+k^)ei(G+k)τjCnk(G)(18)=1VGilgj(|G+k|)Ylm(G+k^)eiGτjCnk(G)

In the first step, one has to find out all the plane-waves within the sphere defined by Eq.(7). Secondly, gj values on the relevant |G+k| are obtained by Spline-interpolation and multiplied with the corresponding real spherical harmonics Ylm. The sum over 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
p~jk|u~nk=|rτj|<rjceik(rτj)fj(|rτ|)Ylm(rτ^)u~nk(r)dr(19)=|rkτj|<rjceikk(rτj)fj(|rkτj|)Ylm(rkτ^j)u~nk(rk)VNfft

First, the grid points within the PAW sphere |rτj|<rjc are identified. Secondly, the values of fj on these grid points are obtained by Spline-interpolation and multiplied with the corresponding real spherical harmonics Ylm. Finally, muliplying the phase term and u~nk(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 ϕj and ϕ~j are defined on the radial logarithmic grid with center at τj (in the first unit cell), i.e.

(20)r|ϕj=ϕj(rτj)(21)r|ϕ~j=ϕ~j(rτj)

The radial logarithmic grid is defined as

rn=r0enH;n=0,1,2,,N1

or equivalently

rn+1=rneH;n=0,1,2,,N1

where r0 is the minimum value of |rτj| and is by definition not zero, N is the number of radial grid points.

Note that rN1 is not necessarily equal to the PAW sphere cutoff radius rjc, as can be seen in Figure 2. However, ϕj and ϕ~j are indential beyond |rτj|=rjc, which makes the difference of the two zero in Eq.(1) beyond rjc.

The AE/PS partial waves are again composed of radial and angular parts,

(22)ϕj(rτj)=Rj(|rτj|)|rτj|Ylm(rτj^)(23)ϕ~j(rτj)=R~j(|rτj|)|rτj|Ylm(rτj^)

IT MUST BE STRESSED that Rj(r) and 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. VASP expands the PS partial waves within the sphere in terms of a linear combination of spherical Bessel functions 5 6

    (24)R~j(r)r={i=12αijl(qkr),r<rcRj(r)r,r>rc

    where the coefficients αi and qi are determined so that the PS partial waves are two times continuously differentiable at rc.

    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 ψnk 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 unk(r)

|unk=|u~nk+eikrj[|ϕj(rτj)|ϕ~j(rτj)]p~j(rτj)|ψ~nk(25)=|u~nk+jeik(rτj)[|ϕj(rτj)|ϕ~j(rτj)]βnkjk

Similar to the evaluation of Eq.(17), Eq.(25) 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 unk. This can be done by introducting a new energy cutoff Eae, which must be larger than the energy cutoff Ec that determines the size of the uniform grid used to represent u~nk. Obviously, the new and old grid size are related by

(26)(nxae,nyae,nzae)=EaeEc(nxc,nyc,nzc)

where denotes a integer ceiling operation.

Secondly, Fourier transform the plane-wave ceofficients Cnk(G) to get u~nk on the new grid.

The rest is to get the real-space representation of the on-site contribution (the summation term in Eq.(25)), which can be evaluated by two approaches.

real-space approach

The on-site term in real space is simply

=jβnkjkeik(rτj)[ϕj(rτj)ϕ~j(rτj)](27)=jβnkjkeik(rτj)[Rj(rτj)R~j(rτj)]|rτj|Ylm(rτj^)

This can be done similar to the the real-space evaluation of βnkjk.

momentum-space approach

Let’s now inspect the plane-wave expansion of the on-site terms

Sj(G)=βnkjkG|eik(rτj)[|ϕj(rτj)|ϕ~j(rτj)]=βnkjkeikτjGG|eikr|GG|[|ϕj(rτj)|ϕ~j(rτj)]=βnkjkeikτjGδG+k,GG|[|ϕj(rτj)|ϕ~j(rτj)]=βnkjkeikτjG+k|[|ϕj(rτj)|ϕ~j(rτj)]=ilVβnkjkeikτj[Gj(|G+k|)G~j(|G+k|)]Ylm(G+k^)ei(G+k)τj(28)=ilVβnkjk[Gj(|G+k|)G~j(|G+k|)]Ylm(G+k^)eiGτj

where G is the spherical Bessel transform of the partical waves

(29)Gj(k)=4π0Rj(r)rjl(kr)r2dr(30)G~j(k)=4π0R~j(r)rjl(kr)r2dr

Note that the plane-waves are now determined by

(31)2|G+k|22me<Eae

The real-space on-site term can be obtained by Fourier transform of jSj(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.(1). 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 Eae=25×Ec.

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

    Figure 3. 2D plot of CO2 HOMO orbital at x = 0 plane. The solid circles indicate the positions of the carbon and oxygen atoms while the dashed circles show the corresponding PAW sphere. The C-O bond length is 1.178A˚ while the PAW cutoff radius for C and O are 0.809A˚ and 0.822A˚, respectively. The related files used to generate this figure can be found in examples of the VaspBandUnfolding package.
  • Wavefunctions ψ(r,z) along the bond direction z at different vertical distances r.

    Figure 4. Comparison of all-electron (blue solid line) and pseudo-wavefunction (red dashed line) of the CO2 HOMO orbital, where z is along the CO2 bond direction and r is the direction perpendicular to the CO2 bond. The circles indicate the positions of the carbon and oxygen atoms. The related files used to generate this figure can be found in examples of the VaspBandUnfolding package.

Reference

This post is licensed under CC BY 4.0 by the author.