8.12. Projector Augmented-Wave Method (PAW)

We can use the Density Functional Theory (DFT) to reduce the many body problem to solve a single particle Schrödinger equation:

H\ket{\psi_n} = \epsilon_n\ket{\psi_n}

The wavefunctions contain cusps (and are oscillatory around each nucleus), also one needs to solve this for all core states.

Next step is to use the known behavior around each atom and take advantage of the known physics somehow. There are two general approaches, either one can incorporate the known physic in the basis (for example the partition of unity in the finite element method), or into the equations. PAW method uses the latter approach.

8.12.1. Projectors, Augmentation Spheres and Smooth Wavefunctions

We introduce smooth wavefunctions (we’ll use “~” for smooth functions) by a linear transformation operator T:

\ket{\psi_n} = T \ket{\tilde\psi_n}

We construct augmentation spheres |{\bf r} - {\bf R}^a| < r_c^a around each atom a (one can imagine a muffin-tin), where r_c^a is a cut-off radius, a is the atom index, {\bf R}^a is the atom position.

We write T as:

T = \one + \sum_a T^a

where T^a only acts in the augmentation sphere. We choose a complete basis \ket{\phi_i^a} (also called partial waves) inside the sphere. The smooth partial waves can be obtained using the T operator:

    = T\ket{\tilde\phi_i^a}
    = \left(\one+\sum_{a'} T^{a'}\right) \ket{\tilde\phi_i^a}

    = (\one+T^a) \ket{\tilde\phi_i^a}

Because T^a only acts in the sphere, it follows that

\ket{\phi_i^a} = \ket{\tilde\phi_i^a}\quad\quad
    \mbox{for $r>r_c^a$}

outside the sphere (i.e. \braket{{\bf r} | \phi_i^a} = \braket{{\bf r} | \tilde\phi_i^a} for r>r_c^a). We can now expand the smooth wavefunctions using the partial waves basis:

(\ket{\tilde\psi_n} = \sum_i P_{ni}^a\ket{\tilde \phi_i^a}

inside the augmentation sphere. Multiplying both sides by T:

(\ket{\tilde\psi_n} = T\sum_i P_{ni}^a\ket{\tilde \phi_i^a}

T\ket{\tilde\psi_n} = \sum_i P_{ni}^a\, T\ket{\tilde \phi_i^a}

\ket{\psi_n} = \sum_i P_{ni}^a \ket{\phi_i^a}

So both smooth and non-smooth wavefunctions have the same expansion coefficients P_{ni}^a. We choose smooth projector functions \ket{\tilde p_i^a} satisfying the following orthogonality and completeness relations inside the augmentation spheres (no restrictions are imposed outside the spheres, so we just define \braket{{\bf r} | \tilde p_i^a} = 0):

(\braket{\tilde p_i^a | \tilde \phi_j^a} = \delta_{ij}

\sum_i \ket{\tilde \phi_i^a}\bra{\tilde p_i^a}  = \one

then multiplying ( by \bra{\tilde p_i^a} and using (

\braket{\tilde p_i^a | \tilde \psi_n }
    = \sum_j P_{nj}^a\braket{\tilde p_i^a | \tilde \phi_j^a}
    = \sum_j P_{nj}^a\delta_{ij}
    = P_{ni}^a

we can rewrite ( and (

(\ket{\tilde\psi_n} = \sum_i \braket{\tilde p_i^a | \tilde \psi_n }
    \ket{\tilde \phi_i^a}

\ket{\psi_n} = \sum_i \braket{\tilde p_i^a | \tilde \psi_n }

Let’s write T^a using the projectors:

    = T^a \one
    = T^a \sum_i \ket{\tilde \phi_i^a}\bra{\tilde p_i^a}
    = \sum_i (T^a\ket{\tilde \phi_i^a})\bra{\tilde p_i^a}
    = \sum_i (\ket{\phi_i^a} - \ket{\tilde \phi_i^a})\bra{\tilde p_i^a}

Note that the right hand side is zero outside the augmentation sphere. Thus

    = \one + \sum_a T^a
    = \one + \sum_a \sum_i (\ket{\phi_i^a} - \ket{\tilde \phi_i^a})\bra{\tilde p_i^a}

In other words, the transformation operator T is completely defined using the smooth and non-smooth partial waves and the projector functions. In terms of the wavefunction:

\ket{\psi_n} = T\ket{\tilde\psi_n}
    =\ket{\tilde\psi_n} + \sum_a \sum_i
        (\ket{\phi_i^a} - \ket{\tilde \phi_i^a})
        \braket{\tilde p_i^a | \tilde\psi_n} =

    =\ket{\tilde\psi_n} + \sum_a\left(
        \sum_i \ket{\phi_i^a}\braket{\tilde p_i^a | \tilde\psi_n}
        - \sum_i\ket{\tilde \phi_i^a}\braket{\tilde p_i^a | \tilde\psi_n}

In words, the wavefunction can be decomposed as the sum of the smooth wavefunction and sum over atoms (centers), at each atom we have “1-center all electron” minus “1-center pseudo”.

The projection functions can always be written as

\bra{\tilde p_i^a}
    = \sum_j\left\{\braket{f_k^a | \tilde\phi_l^a}\right\}_{ij}^{-1}

where \ket{f_k^a} is any set of linearly independent functions.

Note: the n above means all states of interest — either all states, or only the valence states.

8.12.2. Frozen Core Approximation

One can either calculate all electrons in the eigenproblem, or only calculate the valence electrons and treat the core states separately. The simplest option is to introduce a frozen core approximation, where

\ket{\psi_n} = \ket{\phi_\alpha^{a,\mbox{core}}}

for all core states n, here n runs over (a, \alpha), where a is the atom index and \alpha are the core states of an atom. This approximation can also be relaxed in various ways.

8.12.3. Expectation Values of Local Operators

In the frozen core approximation:

\braket{O} = \sum_n^{\mbox{val}} f_n \braket{\psi_n | O | \psi_n}
    + \sum_a \sum_\alpha^{\mbox{core}}
    \braket{\phi_\alpha^{a,\mbox{core}} | O | \phi_\alpha^{a,\mbox{core}}}
    = \cdots =

  =\sum_n^{\mbox{val}} f_n \braket{\tilde \psi_n | O | \tilde \psi_n}
  +\sum_a \sum_{i, j} \left(
    \braket{\phi_i^a | O | \phi_j^a} - \braket{\tilde \phi_i^a | O | \tilde
    \phi_j^a}\right) D_{ij}^a
    \sum_a \sum_\alpha^{\mbox{core}}
    \braket{\phi_\alpha^{a,\mbox{core}} | O | \phi_\alpha^{a,\mbox{core}}}

where the tensor D_{ij}^a is:

D_{ij}^a = \sum_n f_n \braket{\tilde \psi_n|\tilde p_i^a}
    \braket{\tilde p_j^a|\tilde\psi_n}


n({\bf r})
    = \sum_n f_n |\psi_n({\bf r})|^2
    = \sum_n f_n \braket{\psi_n | {\bf r} }\braket{{\bf r} | \psi_n}
    = \Big< \ket{\bf r}\bra{\bf r} \Big> =

  =\sum_n^{\mbox{val}} f_n |\tilde\psi_n({\bf r})|^2
  +\sum_a \sum_{i, j} \left(
    \phi_i^a \phi_j^a - \tilde \phi_i^a \tilde \phi_j^a\right) D_{ij}^a
    \sum_a \sum_\alpha^{\mbox{core}}

The functions \phi_\alpha^{a,\mbox{core}} are not strictly localized withing the augmentation sphere.

8.12.4. Kohn Sham Equations

We multiply the original equations by T^\dag from the left and introduce the smooth wavefunctions:

H\ket{\psi_n} = \epsilon_n\ket{\psi_n}

T^\dag H\ket{\psi_n} = \epsilon_n T^\dag\ket{\psi_n}

T^\dag H T\ket{\tilde\psi_n} = \epsilon_n T^\dag T \ket{\tilde\psi_n}

The orthogonality of wavefunctions is:

\braket{\psi_n | \psi_m} = \delta_{nm}

\braket{\tilde \psi_n | T^\dag T | \tilde \psi_m} = \delta_{nm}

The overlap operator T^\dag T can be written as:

T^\dag T = \cdots = \one + \sum_a \sum_{i,j}
    \ket{\tilde p_i^a} Q_{ij} \bra{\tilde p_j^a}


Q_{ij} = \braket{\phi_i^a | \phi_j^a}
        -\braket{\tilde \phi_i^a | \tilde\phi_j^a}

The transformed Hamiltonian is

H =-\half\nabla^2 + V_H({\bf r}) + V_{xc}({\bf r}) + v({\bf r})

T^\dag H T = \cdots =
    -\half\nabla^2 + V_H(\tilde n) + V_{xc}(\tilde n) +
        \sum_a \sum_{ij}\ket{\tilde p_i^a} H_{ij}^a \bra{\tilde p_j^a}


H_{ij}^a = \braket{\phi_i^a |-\half \nabla^2 + v_{\mbox{eff}} | \phi_j^a}
        -\braket{\tilde \phi_i^a |-\half \nabla^2 + \tilde v_{\mbox{eff}} | \tilde\phi_j^a}