Projected Augmented Wave (PAW) Method: extended resolution of unity method

Abstract

The Projected Augmented Wave (PAW) Method is based on generating a transform between the pseudo wavefunction and all electron wavefunction. For the accuracy of the method, it is important that the local part of the transform (inside each atomic sphere $\mathbf{R}$) be over a complete basis set (with deviations from completeness leading to corrections to the total energy not computed within current implementations of PAW). Here we show how to make this basis much closer to complete without significant additional computational work without modifying the transformation in any way by extending the resolution of unity used for the transform to include more wavefunctions and having them transform via the identity.

I Introduction

It is well know that all electron methods are more technically difficult then pseudopotential methods [3, 8, 7, 13, 12, 1, 11, 9, 10, 2, 14, 6, 5, 4]. This is due to the complex nature of the basis wavefunctions for all electron methods which in particular generically involve augmentation. As such the basis wavefunctions are not of plane wave form, but are augmented inside the Muffin Tin sphere, and therefore require additional calculations for a variety of applications of Density Functional Theory (DFT) methods. In particular:

• The calculation of eigenvalues and eigenvectors of the Khon Sham (KS) Hamiltonian requires FLAPW (Full Potential Linearized Augmented Plane Waves [22, 16, 18, 21, 15, 19, 20, 17]) which requires a relatively complex setup. On the other hand full potential plane wave methods have simple expressions for KS Hamiltonian matrix elements and produce the identity matrix for wavefunction overlaps.

• The Coulomb term (Hartree Energy) is more complex is augmented systems it requires the Weinert method for efficient calculation and convergence without shape approximations [23], where as for pseudopotentials Ewald sums suffice for most applications [24].

• Atomic force calculations in all electron methods are more complex as they require the computation of the Pulay contribution to the forces - which arise due to changes of span of basis set with atomic positions (the Hellman-Feynman theorem does not apply directly in this case [25]).

• Substitutional alloys are difficult to study for all electron methods (other then through cumbersome methods such as Korringa Khon Rostoker (KKR) methods [27, 26]), where as for pseudopotentials there are many viable virtual crystal approximation methods [28].

• The Car Parrinello molecular dynamics simulations [29, 24] method is significantly easier with pseudopotentials then with augmented waves such as Linearize Augmented plane Waves (LAPW) and very limited work has been done with all electron methods in that direction [24, 1].

As such despite many advances in the accuracy of all electron methods [32, 31, 30, 11, 12, 10, 9, 7, 8] it is extremely worthwhile to study pseudopotential methods (atleast currently) and they are currently the leading methods for electronic structure calculations within DFT. Modern pseudopotential methods may be divided into three different categories:

• Norm conserving pseudopotentials [2, 24]

• Ultrasoft psudopotentials [33]

• Projected Augmented Waves (PAW) [34]

We now describe these three methods in turn. Norm conserving pseudopotentials are based on the idea that a perfectly spherical scattering potential can be represented by a scattering phase shift - in each partial wave and at each energy - as a boundary condition at the edge of the sphere. Any potential producing these scattering phase shifts is sufficient for calculations of eigenvalues of the KS Hamiltonian. Therefore an equivalent problem where a much smoother potential with the same scattering phase shift (in the middle of the valence band) as the atomic potential, in each different angular momentum channel, can be used in the KS Hamiltonian. These smoother potentials may be used to solve for the eigenvalues of the KS problem leading to lower cutoffs needed for electronic structure calculations. The norm conserving condition insures both accurate Madelung energy and that the derivative of the scattering phase shift with respect to energy is the same both for the atomic potential and the pseudopotential leading to improved accuracy as a function of energy. On the other hand, ultrasoft potentials use the many state generalizations of the Kleinman Bylander transform [2] in order to study scattering simultaneously at many energies thereby getting rid of the norm conserving condition (needed for accuracy without this step) and restore correct Madelung energies with a a fictitious overlap matrix. PAW introduces a transformation between the soft (low cutoff wavefunctions) and the all electron wavefunctions. PAW is sometimes considered the most versatile because at its core its an efficient method to map a pseudo-wave function problem onto an all electron problem leading to great accuracy of computations. One of its limitations is that it assumes some completeness relations for pseudo wavefunctions used in the transform, a problem we here partially remedy. Indeed we write the all electron wavefunction is given through the relationship

$$\psi=\mathcal{T}\tilde{\psi}$$ (1)

where $\tilde{\psi}$ is the a pseudo wavefunction. We now write:

$$\mathcal{T}=\mathbb{I}+\sum_{\mathbf{R}}S_{\mathbf{R}}$$ (2)

Here $\mathbf{R}$ are the sites of the atomic nuclei. Where

$$S_{\mathbf{R}}=\sum_{i\in\mathbf{R}}\left[\left|\phi_{i}\right\rangle-\left|\tilde{\phi}_{i}\right\rangle\right]\left\langle\tilde{p}_{i}\right|$$ (3)

Where $\left|\phi_{i}\right\rangle$ are from atomic calculations representing the exact atomic states while $\left|\tilde{\phi}_{i}\right\rangle$ are smooth functions. Where

$$\left\langle\tilde{p}_{i}\mid\tilde{\phi}_{j}\right\rangle=\delta_{ij}$$ (4)

There are many options for the wavefunctions $\left|\tilde{p}_{i}\right\rangle$, for example:

$$\left|\tilde{p}_{i}\right\rangle=\sum_{j\in\mathbf{R}}\left|\tilde{\phi}_{j}\right\rangle\left[M\right]_{ji}^{-1}$$ (5)

Where

$$M_{ij}=\left\langle\tilde{\phi}_{j}\mid\tilde{\phi}_{i}\right\rangle$$ (6)

This means that

$$\left[\sum_{i\in\mathbf{R}}\left|\tilde{\phi}_{i}\right\rangle\left\langle\tilde{p}_{i}\right|\right]\left|\tilde{\phi}_{j}\right\rangle=\left|\tilde{\phi}_{j}\right\rangle$$ (7)

so that

$$\sum_{i\in\mathbf{R}}\left|\tilde{\phi}_{i}\right\rangle\left\langle\tilde{p}_{i}\right|=\mathbb{I}_{\mathbf{R}}$$ (8)

provided that the $\left|\tilde{\phi}_{i}\right\rangle$ form a complete basis set. However we have that in reality:

$$\sum_{i\in\mathbf{R}}\left|\tilde{\phi}_{i}\right\rangle\left\langle\tilde{p}_{i}\right|=\left\{\begin{array}[]{cc}\mathbb{I}_{\mathbf{R}}&\in Span\left\{\left|\tilde{\phi}_{i}\right\rangle\right\}\\ 0&otherwise\end{array}\right.$$ (9)

This relation, assumed completeness Eq. (8), is used to prove the identity that for a local operator:

$$\left\langle\tilde{\psi}\right|\hat{A}\left|\tilde{\psi}\right\rangle_{\mathbf{R}}=\left\langle\tilde{\psi}_{\mathbf{R}}\right|\hat{A}\left|\tilde{\psi}_{\mathbf{R}}\right\rangle$$ (10)

Where $\left\langle\right\rangle_{\mathbf{R}}$ means all integrations are done over the atomic sphere centered about $\mathbf{R}$. Where:

$$\left|\psi_{\mathbf{R}}\right\rangle=\left[\sum_{i\in\mathbf{R}}\left|\phi_{i}\right\rangle\left\langle\tilde{p}_{i}\right|\right]\left|\tilde{\psi}\right\rangle$$ (11)

These relations are used to compute the form of the density of electrons the correlation and exchange and the Hartree (Coulomb piece) of the Hamiltonian [34, 2, 24]. This is essential to current PAW methodology.

The main problem is that $\left|\tilde{\phi}_{i}\right\rangle$ do not form a complete basis set because it is of the same cardinality as $\left|\phi_{i}\right\rangle$ and we must obtain $\left|\phi_{i}\right\rangle$ from atomic calculations where typically only two or three wavefunctions, per angular momentum channel, are available with current data sets. Here we propose to circumvent this bottleneck by introducing:

$$\psi=\mathcal{T}_{M}\tilde{\psi}$$ (12)

where $\tilde{\psi}$ is the a pseudo wavefunction and $M$ stands for multiple states. We now write:

$$\mathcal{T}_{M}=\mathbb{I}+\sum_{\mathbf{R}}S_{\mathbf{R}}^{M}$$ (13)

Where

$$S_{\mathbf{R}}^{M}=\sum_{i,\alpha\in\mathbf{R}}\left[\left|\phi_{i\alpha}\right\rangle-\left|\tilde{\phi}_{i\alpha}\right\rangle\right]\left\langle\tilde{p}_{i\alpha}\right|$$ (14)

Where

$$\left|\phi_{i\alpha}\right\rangle=\left\{\begin{array}[]{cc}\left|\phi_{i}\right\rangle&\alpha=1\\ \left|\tilde{\phi}_{i\alpha}\right\rangle&\alpha>1\end{array}\right.$$ (15)

Here $\alpha$ runs over a list of on the order of ten additional pseudo wavefunctions per an all electron wavefunction. Where we have used the idea that it is much easier to generate smooth (low cutoff functions say a combination of Bessel functions with different energies times spherical harmonics with different angular momentum parts) then to do atomic data set calculations. As such multiple smooth wavefunctions can easily be mapped onto the same all electron atomic calculation wavefunction. Where we have that:

$$\left\langle\tilde{p}_{i\alpha}\mid\tilde{\phi}_{j\beta}\right\rangle=\delta_{ij}\delta_{\alpha\beta}$$ (16)

There are many options for the wavefunctions $\left|\tilde{p}_{i\alpha}\right\rangle$ for example:

$$\left|\tilde{p}_{i\alpha}\right\rangle=\sum_{j\in\mathbf{R}}\left|\tilde{\phi}_{j\beta}\right\rangle\left[M\right]_{j\beta i\alpha}^{-1}$$ (17)

Where

$$M_{i\alpha j\beta}=\left\langle\tilde{\phi}_{i\alpha}\mid\tilde{\phi}_{j\beta}\right\rangle$$ (18)

We notice that

$$S_{\mathbf{R}}^{M}=S_{\mathbf{R}}$$ (19)

Except now:

$$\sum_{i\in\mathbf{R}}\left|\tilde{\phi}_{i\alpha}\right\rangle\left\langle\tilde{p}_{i\alpha}\right|=\left\{\begin{array}[]{cc}\mathbb{I}_{\mathbf{R}}&\in Span\left\{\left|\tilde{\phi}_{i\alpha}\right\rangle\right\}\\ 0&otherwise\end{array}\right.$$ (20)

So is much closer to the resolution of identity, as such we have inserted the identity transform in a clever way improving the correctness of the resolution of identity.

II Some technical details

II.1 Main modifications to local operators

We now write:

$$\left|\psi_{\mathbf{R}}^{M}\right\rangle =\left[\sum_{i\alpha\in\mathbf{R}}\left|\phi_{i\alpha}\right\rangle\left\langle\tilde{p}_{i\alpha}\right|\right]\left|\tilde{\psi}\right\rangle$$

$$\left|\tilde{\psi}_{\mathbf{R}}^{M}\right\rangle =\left[\sum_{i\alpha\in\mathbf{R}}\left|\phi_{i\alpha}\right\rangle\left\langle\tilde{p}_{i\alpha}\right|\right]\left|\tilde{\psi}\right\rangle$$ (21)

Now for a local operator $\hat{A}$ we write:

$$\left\langle\mathcal{T}_{M}^{\dagger}\hat{A}\mathcal{T}_{M}\right\rangle$$

$$=\left[\left\langle\tilde{\psi}\right|+\sum_{\mathbf{R}}\left[\left\langle\psi_{\mathbf{R}}^{M}\right|-\left\langle\tilde{\psi}_{\mathbf{R}}^{M}\right|\right]\right]\times$$

$$\times\hat{A}\left[\left|\tilde{\psi}\right\rangle+\sum_{\mathbf{R}}\left[\left|\psi_{\mathbf{R}}^{M}\right\rangle-\left|\tilde{\psi}_{\mathbf{R}}^{M}\right\rangle\right]\right]$$

$$=\left[\left\langle\tilde{\psi}\right|-\sum_{\mathbf{R}}\left\langle\tilde{\psi}_{\mathbf{R}}^{M}\right|\right]\hat{A}\left[\left|\tilde{\psi}\right\rangle-\sum_{\mathbf{R}^{\prime}}\left|\tilde{\psi}_{\mathbf{R}^{\prime}}^{M}\right\rangle\right]+$$

$$+\sum_{\mathbf{R}}\left\langle\psi_{\mathbf{R}}^{M}\right|\hat{A}\left|\psi_{\mathbf{R}}^{M}\right\rangle$$

$$=\left\langle\tilde{\psi}\right|\hat{A}\left|\tilde{\psi}\right\rangle-\sum_{\mathbf{R}}\left\langle\tilde{\psi}_{\mathbf{R}}^{M}\right|\hat{A}\left|\tilde{\psi}_{\mathbf{R}}^{M}\right\rangle+\sum_{\mathbf{R}}\left\langle\psi_{\mathbf{R}}^{M}\right|\hat{A}\left|\psi_{\mathbf{R}}^{M}\right\rangle$$ (22)

It is now also useful to introduce:

$$D_{i\alpha j\beta}^{\mathbf{R}}=\sum_{n,\mathbf{k}}f\left(\varepsilon_{n}\left(\mathbf{k}\right)\right)\left\langle\tilde{\psi}_{n}\left(\mathbf{k}\right)\right|\left[\left|\tilde{p}_{i\alpha}\right\rangle\left\langle\tilde{p}_{j\beta}\right|\right]\left|\tilde{\psi}_{n}\left(\mathbf{k}\right)\right\rangle$$ (23)

This means that:

$$\left\langle\mathcal{T}_{M}^{\dagger}\hat{A}\mathcal{T}_{M}\right\rangle$$

$$=\sum_{n,\mathbf{k}}f\left(\varepsilon_{n}\left(\mathbf{k}\right)\right)\left\langle\tilde{\psi}_{n}\left(\mathbf{k}\right)\right|\hat{A}\left|\tilde{\psi}_{n}\left(\mathbf{k}\right)\right\rangle+$$

$$+\sum_{\mathbf{R},i,j}D_{i\alpha j\beta}^{\mathbf{R}}\left\langle\phi_{i\alpha}\right|\hat{A}\left|\phi_{j\beta}\right\rangle-\sum_{\mathbf{R},i\alpha,j\beta}D_{i\alpha j\beta}^{\mathbf{R}}\left\langle\tilde{\phi}_{i\alpha}\right|\hat{A}\left|\tilde{\phi}_{j\beta}\right\rangle$$

$$+\sum_{core\in\mathbf{R}}\left\langle\phi_{c,i}\right|\hat{A}\left|\phi_{c,i}\right\rangle$$ (24)

Here we have introduced core states from the atomic calculations. Notice that this is a different expression then the one with regular PAW in particular there are cross terms now between $\alpha=1$ and $\alpha>1$. In particular works for

$$\hat{A}=\delta\left(\mathbf{r}-\mathbf{r}^{\prime}\right),\>\hat{A}=\frac{-\nabla^{2}}{2m}$$ (25)

which are important for ground state energy determination. That is for the density and kinetic energy terms. In particular:

$$\mathcal{T}_{M}^{\dagger}\hat{H}_{KS}\mathcal{T}_{M}\left|\tilde{\psi}_{n}\right\rangle=\varepsilon_{n}\mathcal{T}_{M}^{\dagger}\mathcal{T}_{M}\left|\tilde{\psi}_{n}\right\rangle$$ (26)

Indeed

$$S_{M} =\mathcal{T}_{M}^{\dagger}\mathcal{T}_{M}=1+\sum_{\mathbf{R}}S_{\mathbf{R}}^{M}+\sum_{\mathbf{R}}S_{\mathbf{R}}^{M\dagger}+\sum_{\mathbf{R}}S_{\mathbf{R}}^{M\dagger}S_{\mathbf{R}}^{M}$$

$$=1+\left[\sum_{i\alpha\in\mathbf{R}}\left|\tilde{p}_{i\alpha}\right\rangle\left\langle\tilde{\phi}_{i\alpha}\right|\right]\sum_{i,\alpha\in\mathbf{R}}\left[\left|\phi_{j\beta}\right\rangle-\left|\tilde{\phi}_{j\beta}\right\rangle\right]\left\langle\tilde{p}_{j\beta}\right|+$$

$$+\sum_{i,\alpha\in\mathbf{R}}\left|\tilde{p}_{i\alpha}\right\rangle\left[\left\langle\phi_{i\alpha}\right|-\left\langle\tilde{\phi}_{i\alpha}\right|\right]\left[\sum_{j\beta\in\mathbf{R}}\left|\tilde{\phi}_{j\beta}\right\rangle\left\langle\tilde{p}_{j\beta}\right|\right]+$$

$$+\sum_{i,\alpha\in\mathbf{R}}\left|\tilde{p}_{i\alpha}\right\rangle\left[\left\langle\phi_{i\alpha}\right|-\left\langle\tilde{\phi}_{i\alpha}\right|\right]\sum_{j\beta\in\mathbf{R}}\left[\left|\phi_{j\beta}\right\rangle-\left|\tilde{\phi}_{j\beta}\right\rangle\right]\left\langle\tilde{p}_{j\beta}\right|$$

$$=1+\sum_{i,\alpha\in\mathbf{R}}\sum_{j\beta\in\mathbf{R}}\left|\tilde{p}_{i\alpha}\right\rangle\left[\left\langle\phi_{i\alpha}\mid\phi_{j\beta}\right\rangle-\left\langle\tilde{\phi}_{i\alpha}\mid\tilde{\phi}_{j\beta}\right\rangle\right]\left\langle\tilde{p}_{j\beta}\right|$$

$$\equiv 1+\sum_{i,\alpha\in\mathbf{R}}\sum_{j\beta\in\mathbf{R}}\left|\tilde{p}_{i\alpha}\right\rangle\Delta_{i\alpha j\beta}\left\langle\tilde{p}_{j\beta}\right|\neq\mathcal{T}^{\dagger}\mathcal{T}=S$$ (27)

The reason for the inequality is that we have used a different (more accurate) resolution of identity to simplify the situation which leads to different results. Similarly for

$$\mathcal{T}_{M}^{\dagger}\hat{H}_{KS}\mathcal{T}_{M}=\hat{H}_{KS}+$$

$$\sum_{i,\alpha\in\mathbf{R}}\sum_{j\beta\in\mathbf{R}}\left|\tilde{p}_{i\alpha}\right\rangle\left[\left\langle\phi_{i\alpha}\mid\hat{H}_{KS}\mid\phi_{j\beta}\right\rangle-\left\langle\tilde{\phi}_{i\alpha}\mid\hat{H}_{KS}\mid\tilde{\phi}_{j\beta}\right\rangle\right]\left\langle\tilde{p}_{j\beta}\right|$$

$$\equiv\hat{H}_{KS}+\sum_{i,\alpha\in\mathbf{R}}\sum_{j\beta\in\mathbf{R}}\left|\tilde{p}_{i\alpha}\right\rangle h_{i\alpha j\beta}\left\langle\tilde{p}_{j\beta}\right|\neq\mathcal{T}^{\dagger}\hat{H}_{KS}\mathcal{T}$$ (28)

Again the resolution of identity changed the final result.

II.2 Correlation and exchange and Hartree piece

We write for local correlation and exchange functionals such as Local Density Approximation (LDA) or Generalized Gradient Approximation (GGA):

$$E_{XC}=E_{XC}\left(\tilde{n}\right)+\sum_{\mathbf{R}}E_{XC}\left(n_{\mathbf{R}}\right)-\sum_{\mathbf{R}}E_{XC}\left(\tilde{n}_{\mathbf{R}}\right)$$ (29)

This is based on the idea that

$$\tilde{n}=\tilde{n}_{\mathbf{R}}\>\left(in\>the\>sphere\>\mathbf{R}\right)$$ (30)

However for this we see that Eq. (20) is much better then Eq. (9) with explicit changes in the form of $\tilde{n}_{\mathbf{R}}$ seen from Eq. (24). Furthermore we now introduce the convenient notation:

$$\left(f\mid g\right)=\int d^{3}\mathbf{r}_{1}\int d^{3}\mathbf{r}_{2}\frac{f\left(\mathbf{r}_{1}\right)g\left(\mathbf{r}_{2}\right)}{\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|}$$ (31)

and

$$\left(f\mid f\right)=\left(\left(f\right)\right)$$ (32)

Similarly we have that the Coulomb piece is given by:

$$E_{C}\left[n\right] =U_{H}\left(\tilde{\rho}\right)$$

$$+\frac{1}{2}\sum_{\mathbf{\mathbf{R}}}\left[\left(\left(n_{\mathbf{R}}\right)\right)+2\left(n_{\mathbf{R}}\mid Z_{\mathbf{R}}\right)-\left(\left(\tilde{n}_{\mathbf{R}}+\tilde{Z}_{\mathbf{R}}\right)\right)\right]$$ (33)

where

$$U_{H}\left(\tilde{\rho}\right)=\frac{1}{2}\left(\left(\tilde{n}+\tilde{Z}_{\mathbf{R}}\right)\right)$$ (34)

Where $Z_{\mathbf{R}}$ is the coulomb charge density and $\tilde{Z}_{\mathbf{R}}$ is the pseudized charge density satisfying the relationship:

$$\int d^{3}\mathbf{r}\left|\mathbf{r}-\mathbf{R}\right|^{l}Y_{lm}\left(\widehat{\mathbf{r}-\mathbf{R}}\right)\times$$

$$\times\left(n_{\mathbf{R}}\left(\mathbf{r}\right)-\tilde{n}_{\mathbf{R}}\left(\mathbf{r}\right)+Z_{\mathbf{R}}\left(\mathbf{r}\right)-\tilde{Z}_{\mathbf{R}}\left(\mathbf{r}\right)\right)$$

$$=0$$ (35)

Where we have repeatedly used Eq. (30) so Eq. (20) is much better then Eq. (9).

III Conclusion

In this work we have introduced a new version of the PAW transformation (given in Eqs. (12), (13), (14) and (15)) which allows for an extended resolution of identity to map between the pseudized problem and the all electron problem. This greatly improves the accuracy with which mapping reproduces the key PAW equation (given in Eq. (10)) as well as the various applications of that identity to exchange and correlation energies, Coulomb energies and the KS secular equation. This improves the accuracy of PAW at very limited computational cost. In the future it would be of interest to include strong correlations into the calculations in the form of LDA+U or LDA+DMFT [35] with this new PAW method.