Requests for technical support from the VASP group should be posted in the VASP-forum.

RPA/ACFDT: Correlation energy in the Random Phase Approximation

ACFDT stands for the adiabatic connection fluctuation dissipation theorem and is an alternative way to derive the energy expression for the correlation energy in the random phase approximation. In the following the diagrammatic description is presented. For the ACFDT formulation the reader is referred to the literature.[1]

Diagrammatic approach to the correlation energy

The correlation energy ${\displaystyle E_{c}}$ is defined as the missing piece of the Hartree-Fock energy ${\displaystyle E_{x}}$ to the total energy, that is ${\displaystyle E_{tot}=E_{x}+E_{c}}$. The exact form of ${\displaystyle E_{c}}$ is unknown and can be calculated only approximately for a realistic system. The Random Phase Approximation (RPA) is such an approximation that provides access to ${\displaystyle E_{c}}$. The RPA was first studied by Bohm and Pines for the homogeneous electron gas and was later recognized by Gell-Mann and Brueckner as an approximation of ${\displaystyle E_{c}}$ that can be expressed in the same language as Feynman used a few years earlier to describe the positron. [2][3][4]

Feynman's diagrammatic approach is based on quantum field theory (QFT), which in turn is based on the Gell-Mann and Low theorem. This theorem states that the eigenstate of an interacting Hamiltonian can be expressed in terms of the eigenstates of the non-interacting one.[5] For this reason, each diagrammatic calculation, like the RPA or GW, requires the solution of the non-interacting Hamiltonian ${\displaystyle H_{0}}$ of the system, like for instance the Hartree-Fock energies and orbitals or the solutions of the Kohn-Sham Hamiltonian ${\displaystyle \epsilon _{n{\bf {k}}},\phi _{n{\bf {k}}}}$.

QFT is commonly formulated in the Dirac (also known as interaction) picture, where dynamics described by the interaction part ${\displaystyle {\hat {V}}}$ of the fully interacting Hamiltonian ${\displaystyle {\hat {H}}={\hat {H}}_{0}+{\hat {V}}}$ are singled out via time-dependent operators like ${\displaystyle {\hat {V}}(t)=e^{i{\hat {H}}_{0}t}{\hat {V}}e^{-i{\hat {H}}_{0}t}}$. These operators act on states like the non-interacting groundstate of the system ${\displaystyle |\Psi _{0}\rangle }$ causing fluctuations at a specific point in time. The main idea of QFT is to understand observations, which can be measured by an observer, as a collective phenomenon of all possible fluctuations.[6]

Thereby, fluctuations are understood as the creation of virtual electrons (and holes) that interact with each other and are annihilated after some time. Formally this is achieved by introducing creation ${\displaystyle {\hat {\psi }}^{\dagger }({\bf {r}},t)}$ and annihilation operators ${\displaystyle {\hat {\psi }}({\bf {r}},t)}$ that satisfy following relations

${\displaystyle {\hat {\psi }}({\bf {r}},t)|\Psi _{0}\rangle =0=\langle \Psi _{0}|{\hat {\psi }}({\bf {r}},t)}$

${\displaystyle \lbrace \psi ({\bf {r}},t),\psi ^{\dagger }({\bf {r}},t)\rbrace =\psi ^{\dagger }({\bf {r}},t)\psi ^{\dagger }({\bf {r}},t)+\psi ^{\dagger }({\bf {r}},t),\psi ({\bf {r}},t)=i\delta ({\bf {r}}-{\bf {r}}')}$

${\displaystyle \lbrace \psi ^{\dagger }({\bf {r}},t),\psi ^{\dagger }({\bf {r}},t)\rbrace =0=\lbrace \psi ({\bf {r}},t),\psi ({\bf {r}},t)\rbrace .}$

The first relation defines the non-interacting groundstate ${\displaystyle |\Psi _{0}\rangle }$ as the Fermi vacuum (the groundstate in the absence of any fluctuations), while the second and third anti-commutator relations are a consequence of the Pauli principle. In fact, all operators that describe measurable quantities of a system of interacting electrons can be represented in terms of ${\displaystyle \psi ^{\dagger }({\bf {r}},t)}$ and ${\displaystyle \psi ({\bf {r}},t)}$ alone; additional objects are not necessary.

However, the time-ordering operator

${\displaystyle {\hat {T}}{\hat {A}}(t){\hat {B}}(t')=\Theta (t-t'){\hat {A}}(t){\hat {B}}(t')-\Theta (t'-t){\hat {B}}(t'){\hat {A}}(t),}$

where ${\displaystyle \Theta (t)}$ is the unit step function, and the time-evolution operator

${\displaystyle {\hat {S}}(t,t_{0})={\hat {T}}e^{-i\int _{t_{0}}^{t}{\hat {V}}(t'){\rm {d}}t'}}$

are helpful quantities, since they allow to formulate the Gell-Mann and Low theorem as follows.

Gell-Mann and Low theorem

Using adiabatic coupling of the interaction ${\displaystyle {\hat {V}}(t)\to {\hat {V}}_{\eta }(t)=e^{-\eta |t|}{\hat {V}}}$, Gell-Mann and Low proved that the vectors

${\displaystyle {\frac {|\Omega _{\nu }\rangle }{\langle \Omega _{\nu }|\Psi _{\nu }\rangle }}=\lim _{\eta \to 0}{\frac {{\hat {S}}_{\eta }(0,-\infty )|\Psi _{\nu }\rangle }{\langle \Omega _{\nu }|\Psi _{\nu }\rangle }}}$

are the eigenstates of the interacting Hamiltonian.[5]

We follow the common literature and suppress the infinitesimal ${\displaystyle \eta }$ in the following bearing in mind that the limit ${\displaystyle \eta \to 0}$ is performed at the very end of the calculation.[7][8]

Diagrammatic perturbation theory

A consequence of the Gell-Mann and Low theorem, is the following form of the interacting groundstate energy[8]

${\displaystyle E_{tot}=E_{0}=\langle \Omega _{0}|{\hat {H}}|\Omega _{0}\rangle ={\frac {\langle \Psi _{0}|{\hat {S}}(\infty ,-\infty ){\hat {H}}|\Psi _{0}\rangle }{\langle \Psi _{0}|{\hat {S}}(\infty ,-\infty )|\Psi _{0}\rangle }},}$

which can be seen as starting point of diagrammatic perturbation theory. The expression above is used to derive all possible approximations by expanding the time-evolution operator ${\displaystyle {\hat {S}}}$ into a series. The resulting matrix-elements of creation and annihilation operators are evaluated term by term using the canonical anti-commutator relations defined above (Wick's theorem[9]). It follows that all terms in perturbation theory are expressed by only two quantities, the non-interacting Feynman propagator

${\displaystyle G_{0}(1,2)=-i\sum _{n{\bf {k}}}\phi ({\bf {r}}_{2})\phi ^{*}({\bf {r}}_{1})e^{-i(\epsilon _{n{\bf {k}}}-\epsilon _{F})(t_{2}-t_{1})}\left[f_{n{\bf {k}}}\Theta (t_{2}-t_{1})-(1-f_{n{\bf {k}}})\Theta (t_{1}-t_{2})\right],\quad 1=({\bf {r}}_{1},t_{1}),2=({\bf {r}}_{2},t_{2})}$

and the Coulomb interaction

${\displaystyle V(1,2)={\frac {\delta (t_{1}-t_{2})}{|{\bf {r}}_{1}-{\bf {r}}_{2}|}}.}$

Then, each term in the series corresponds to an integral over space-time coordinates ${\displaystyle ({\bf {r}},t)}$.

Feynman diagrams are used to illustrate which terms are considered in the perturbation series. The illustration is usually achieved with so-called Feynman rules that map a specific diagram to an integral (and vice versa). For instance the second order diagram

is also known as the direct Møller-Plessett term and stands for following integral

${\displaystyle E_{\rm {dMP}}^{(2)}=\int {\rm {d}}(1,2,3,4)G_{0}(1,2)G_{0}(2,1)V(1,3)V(2,4)G_{0}(3,4)G_{0}(4,3),\quad {\rm {d}}(1,\cdots ,4)={\rm {d}}{\bf {r}}_{1}{\rm {d}}t_{1}\cdots {\rm {d}}{\bf {r}}_{4}{\rm {d}}t_{4}}$

All Feynman rules can be found in the book of Negele and Orland or elsewhere.[7][8]

The Random Phase Approximation

The RPA is obtained from neglecting all second and higher order terms in the perturbation series of the groundstate energy, except of those which can be expressed soley in terms of the independent particle polarizability

${\displaystyle \chi _{0}(1,2)=-iG_{0}(1,2)G_{0}(2,1)}$

corresponding to the "bubble" diagram

Because of the symmetric time property ${\displaystyle \chi _{0}(t_{2}-t_{1})=\chi _{0}(t_{1}-t_{2})}$, the independent particle polarizability is of bosonic character. Because the RPA neglects all non-bosonic terms in the perturbation series, it corresponds essentially to a "bosonization" of the many-body problem for which the n-th order term can be written analytically as[7]

${\displaystyle E_{\rm {dMP}}^{(n)}={\frac {1}{2n}}\int _{-\infty }^{\infty }{\frac {{\rm {d}}\omega }{2\pi }}{\rm {Tr}}\left[{\tilde {\chi }}_{0}(\omega )\cdot V\right]^{n}.}$

Here, the trace of the matrix product is most effectively done in reciprocal space ${\displaystyle \left[{\tilde {\chi }}_{0}(\omega )\cdot V\right]({\bf {q+G}}_{1},{\bf {q+G}}_{2})=\sum _{\bf {G}}{\tilde {\chi }}_{0}({\bf {q+G}}_{1},{\bf {G}},\omega )V({\bf {q+G}},{\bf {q+G}}_{2})}$ using the Fourier transformed polarizability ${\displaystyle {\tilde {\chi }}_{0}({\bf {q+G}}_{1},{\bf {q+G}}_{2},\omega )}$, the diagonal Coulomb potential ${\displaystyle V({\bf {q+G}}_{1},{\bf {q+G}}_{2})={\frac {\delta _{{\bf {G}}_{1}{\bf {G}}_{2}}}{|{\bf {q+G}}_{1}|}}}$ and the conserved crystal momentum ${\displaystyle {\bf {q}}}$ in the first Brillouin zone.

All bubble terms of order ${\displaystyle n\geq 2}$ can be written in a closed form using the series for the logarithm ${\displaystyle \ln(1-x)+x=-\sum _{n=2}^{\infty }{\frac {x^{n}}{n}}}$ and define the correlation part of the RPA energy

${\displaystyle E_{c}^{\rm {RPA}}=\int {\frac {{\rm {d}}\omega }{2\pi }}{\rm {Tr}}\left\lbrace \ln \left[1-{\tilde {\chi }}_{0}(\omega )\cdot V\right]+{\tilde {\chi }}_{0}(\omega )\cdot V\right\rbrace .}$

There are two first order contributions to the total energy that yield the exact exchange energy ${\displaystyle E_{x}=T+V_{ext}+V_{h}+V_{x}}$, which is usually determined separately.

Computational Complexity

The calculation of the RPA integral requires the determination of the independent particle polarizability matrix ${\displaystyle {\tilde {\chi }}_{\bf {GG'}}^{0}({\bf {q}},\omega _{n})={\tilde {\chi }}_{0}({\bf {q+G}},{\bf {q+G}}',\omega _{n})}$ on each of the ${\displaystyle N_{\bf {q}}}$ sampling points of the first Brillouin zone for ${\displaystyle N_{\omega }}$ frequency points. The number of frequency points is reduced drastically, by performing the integration over the imaginary frequency axis ${\displaystyle \omega \to i\omega }$.[10]

The independent particle polarizability on the imaginary axis can be determined with two alternative methods.

Quartic scaling RPA: Direct calculation

Direct calculation of ${\displaystyle {\tilde {\chi }}^{0}}$ using the formula of Adler and Wiser[11] [12]

${\displaystyle {\tilde {\chi }}_{{\bf {GG}}'}^{0}({\bf {q}},i\omega )=\sum \limits _{{\bf {k}}\in BZ}\sum \limits _{n,n'}{\frac {f_{n{\bf {k}}}(1-f_{n{\bf {k-q}}})}{\epsilon _{n{\bf {k-q}}}-\epsilon _{n{\bf {k}}}-i\omega }}\langle \phi _{n{\bf {k-q}}}|e^{i{\bf {Gr}}}|\phi _{n'{\bf {k}}}\rangle \langle \phi _{n'{\bf {k}}}|e^{-i{\bf {G'r'}}}|\phi _{n'{\bf {k-q}}}\rangle }$

yields an RPA algorithm that has a computational cost of ${\displaystyle N_{\omega }N_{\bf {k}}^{2}N_{\bf {G}}^{4}}$. Because the number of plane waves ${\displaystyle N_{\bf {G}}}$ scales linearly with the system size (number of electrons in the unit cell), the direct calculation of the polarizability is unfavourable for large system sizes, e.g. for more than ~20 atoms in the unit cell.

Cubic scaling RPA: Contraction of imaginary time Green's functions

An alternative way to determine ${\displaystyle {\tilde {\chi }}^{0}}$ is to frist determine imaginary time Green's functions of the form[13]

${\displaystyle G_{0}({\bf {r,r'}},i\tau )=\sum \limits _{{\bf {k}}\in BZ}\sum \limits _{n}\phi _{n{\bf {k}}}({\bf {r}})\phi _{n{\bf {k}}}^{*}({\bf {r'}})e^{-(\epsilon _{n{\bf {k}}}-\epsilon _{F})\tau }\left[\Theta (-\tau )f_{n{\bf {k}}}-\Theta (\tau )(1-f_{n{\bf {k}}})\right]}$

and to perform afterwards a Fourier transformation into reciprocal and imaginary frequency space of

${\displaystyle \chi _{0}({\bf {r,r'}},i\tau )=-G_{0}({\bf {r,r'}},i\tau )G_{0}({\bf {r',r}},-i\tau ).}$

Although more evolved, this approach has the advantage that the computational cost for the determination of ${\displaystyle {\tilde {\chi }}_{0}}$ scales with ${\displaystyle N_{\omega }N_{\bf {k}}N_{\bf {G}}^{3}}$ and is essentially only cubic in system size. The space-time method allows to study relatively large systems with the RPA.[14]

Basis set convergence of RPA-ACFDT calculations

The expression for the ACFDT-RPA correlation energy written in terms of reciprocal lattice vectors reads:

${\displaystyle E_{\rm {c}}^{\rm {RPA}}=\int _{0}^{\infty }{\frac {\mathrm {d} \omega }{2\pi }}\sum _{{\mathbf {q} }\in \mathbf {BZ} }\sum _{\mathbf {G} }\left\{(\mathrm {ln} [1-{\tilde {\chi }}^{0}({\mathbf {q} },\mathrm {i} \omega )V({\mathbf {q} })])_{\mathbf {G,G} }+V_{\mathbf {G,G} }({\mathbf {q} }){\tilde {\chi }}^{0}({\mathbf {q} },{\mathrm {i} }\omega )\right\}}$.

The sum over reciprocal lattice vectors has to be truncated at some ${\displaystyle \mathbf {G} _{\mathrm {max} }}$, determined by ${\displaystyle {\frac {\hbar ^{2}|{\mathbf {G} }+{\mathbf {q} }|^{2}}{2\mathrm {m} _{e}}}}$ < ENCUTGW, which can be set in the INCAR file. The default value is ${\displaystyle {\frac {2}{3}}\times }$ ENCUT, which experience has taught us not to change. For systematic convergence tests, instead increase ENCUT and repeat steps 1 to 4, but be aware that the "maximum number of plane-waves" changes when ENCUT is increased. Note that it is virtually impossible, to converge absolute correlation energies. Rather concentrate on relative energies (e.g. energy differences between two solids, or between a solid and the constituent atoms).

Since correlation energies converge very slowly with respect to ${\displaystyle \mathbf {G} _{\rm {max}}}$, VASP automatically extrapolates to the infinite basis set limit using a linear regression to the equation: [1][15][16]

${\displaystyle E_{\mathrm {c} }({\mathbf {G} })=E_{\mathrm {c} }(\infty )+{\frac {A}{{\mathbf {G} }^{3}}}}$.

Furthermore, the Coulomb kernel is smoothly truncated between ENCUTGWSOFT and ENCUTGW using a simple cosine like window function (Hann window function). The default for ENCUTGWSOFT is 0.8${\displaystyle \times }$ENCUTGW (again we do not recommend to change this default).

The integral over ${\displaystyle \omega }$ is evaluated by means of a highly accurate minimax integration.[10] The number of ${\displaystyle \omega }$ points is determined by the flag NOMEGA, whereas the energy range of transitions is determined by the band gap and the energy difference between the lowest occupied and highest unoccupied one-electron orbital. VASP determines these values automatically (from vasp.5.4.1 on), and the user should only carefully converge with respect to the number of frequency points NOMEGA. A good choice is usually NOMEGA=12, however, for large gap systems one might obtain ${\displaystyle \mu }$eV convergence per atom already using 8 points, whereas for metals up to NOMEGA=24 frequency points are sometimes necessary, in particular, for large unit cells.

Strictly adhere to the steps outlines above. Specifically, be aware that steps two and three require the WAVECAR file generated in step one, whereas step four requires the WAVECAR and WAVEDER file generated in step three (generated by setting LOPTICS=.TRUE.).

Matsubara Formalism: Metallic systems at finite Temperature

The zero-temperature formalism of many-body perturbation theory breaks down for metals (systems with zero energy band-gap) as pointed out by Kohn and Luttinger.[17] This conundrum is lifted by considering diagrammatic perturbation theory at finite temperature ${\displaystyle T>0}$, which may be understood by an analytical continuation of the real-time ${\displaystyle t}$ to the imaginary time axis ${\displaystyle -i\tau }$. Matsubara has shown that this Wick-rotation in time ${\displaystyle t\to -i\tau }$ reveals an intriguing connection to the inverse temperature ${\displaystyle \beta =1/T}$ of the system.[18] More precisely, Matsubara has shown that all terms in perturbation theory at finite temperature can be expressed as integrals of imaginary time quantities (such as the polarizability ${\displaystyle \chi (-i\tau )}$) over the fundamental interval ${\displaystyle -\beta \leq \tau \leq \beta }$.

As a consequence, one decomposes imaginary time quantities into a Fourier series with period ${\displaystyle \beta }$ that determines the spacing of the Fouier modes. For instance the imaginary polarizability can be written as

${\displaystyle \chi (-i\tau )={\frac {1}{\beta }}\sum _{m=-\infty }^{\infty }{\tilde {\chi }}(i\nu _{m})e^{-i\nu _{m}\tau },\quad \nu _{m}={\frac {2m}{\beta }}\pi }$

and the corresponding random phase approximation of the correlation energy at finite temperature becomes a series over (in this case bosonic) Matsubara frequencies

${\displaystyle \Omega _{c}^{\rm {RPA}}={\frac {1}{2}}{\frac {1}{\beta }}\sum _{m=-\infty }^{\infty }{\rm {Tr}}\left\lbrace \ln \left[1-{\tilde {\chi }}(i\nu _{m})V\right]-{\tilde {\chi }}(i\nu _{m})V\right\rbrace ,\quad \nu _{m}={\frac {2m}{\beta }}\pi }$

The Matsubara formalism has the advantage that all contributions to the Green's function and the polarizability are mathematically well-defined, including contributions from states close to the chemical potential ${\displaystyle \epsilon _{n{\bf {k}}}\approx \mu }$, such that Matsubara series converge also for metallic systems.

Although formally convenient, Matsubara series converge poorly with the number of considered terms in practice. The VASP code, therefore, uses a compressed representation of the Fourier modes by employing the Minimax Isometry method.[19] This approach converges exponentially with the number of considered frequency points.