ACFDT/RPA calculations

The ACFDT-RPA groundstate energy ( $E_{\mathrm{RPA}}$ ) is the sum of the ACFDT-RPA correlation energy $E_{c}$ and the Hartree-Fock energy evaluated non self-consistently using DFT orbitals $E_{\mathrm{EXX}}$ :

$E_{\mathrm{RPA}}=E_{\mathrm{c}}+E_{\mathrm{EXX}}$ .

Note that, $E_{\mathrm{EXX}}$ here includes also the Hartree energy, the kinetic energy, as well as the Ewald energy of the ions, whereas often in literature $E_{\mathrm{EXX}}$ refers only to the exact exchange energy evaluated using DFT orbitals.

If ALGO=ACFDT is set in the INCAR file, VASP calculates the correlation energy in the random phase approximation. To this end, VASP calculates first the independent particle response function, using the virtual (unoccupied) states found in the WAVECAR file, and then determines the correlation energy using the plasmon fluctuation equation:

$E_{c} = \frac{1}{2 \pi} \int_{0}^{\infty} \mathrm{Tr} \, \ln(1-\chi( i \omega) v) + \chi( i \omega) v \, d \omega$ .

General Recipe to Calculate ACFDT-RPA Total Energies

In practice, RPA energy calculations need to proceed in four steps (VASP can not yet perform all required steps in a single VASP run).

First step (a standard DFT run): All occupied orbitals (and as usual in VASP, a few unoccupied orbitals) of the DFT-Hamiltonian are calculated:

EDIFF = 1E-8
ISMEAR = 0 ; SIGMA = 0.05

This can be done with your favorite setup, but we recommend to attain very high precision (small EDIFF flag) and to use a small smearing width (SIGMA flag), and to avoid higher order Methfessel-Paxton smearing (see also ISMEAR). We suggest to use PBE orbitals as input for the ACFDT-RPA run, but other choices are possible as well, e.g. LDA or hybrid functionals such as HSE. For hybrid functionals, we suggest to carefully consider the caveats mentioned in reference ^[1], specifically the RPA dielectric matrix yields significantly weak screening for hybrid functionals, which deteriorates RPA results.

Second step: the Hartree Fock energy $E_{\mathrm{EXX}}$ is calculated using the predetermined DFT orbitals:

ALGO  = EIGENVAL ; NELM = 1
LWAVE=.FALSE.                  ! avoid accidental update of WAVECAR
LHFCALC = .TRUE. ; AEXX = 1.0  ! you my set ALDAC = 0.0 but the default is 1-AEXX
ISMEAR = 0 ; SIGMA = 0.05

For insulators and semiconductors with a sizable gap, faster convergence of the Hartree-Fock energy can be obtained by setting HFRCUT=-1, altough this slows down k-point convergence for metals.

Third step: Search for maximum number of plane-waves: in the OUTCAR file of the first step, and run VASP again with the following INCAR file to determine all virtual states by an exact diagonalization of the Hamiltonian (DFT or hybrid, make certain to use the same Hamiltonian as in step 1):

NBANDS = maximum number of plane-waves (times 2 for gamma-only calculations)
ALGO = Exact    ! exact diagonalization
NELM = 1        ! one step suffices since WAVECAR is pre-converged
LOPTICS = .TRUE.
ISMEAR = 0 ; SIGMA = 0.05

For calculations using the gamma-point only version of vasp, NBANDS must be set to twice the maximum number of plane-waves: (found in the OUTCAR file) in step 1. For metals, we recommend to avoid setting LOPTICS=.TRUE., since this slows down k-point convergence.^[2]

Fourth step: Calculate the ACFDT-RPA correlation energy:

NBANDS =  maximum number of plane-waves
ALGO = ACFDT
NOMEGA = 8-24

To reach technical convergence, a number of flags are available to control the evaluation of the ACFDT-RPA correlation energy in the fourth step. The expression for the ACFDT-RPA correlation energy reads:

$\mathrm{E}_{\mathrm{c}}=\int_{0}^{\infty} \frac{\mathrm{d}\omega}{2\pi} \sum_{{\mathbf{q}}\in \mathbf{BZ} }\sum_{{\mathbf{G}}} \left\{(\mathrm{ln}[1-\chi^{KS}({\mathbf{q}},\mathrm{i}\omega)\nu({\mathbf{q}})])_{{\mathbf{G,G}}} +\nu_{{\mathbf{G,G}}}({\mathbf{q}})\chi^{KS}({\mathbf{q}},{\mathrm{i}}\omega) \right\}$ .

The sum over reciprocal lattice vectors has to be truncated at some $\mathbf{G}_{\mathrm{max}}$ , determined by $\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 $\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 $\mathbf{G}_{\mathrm<max>}$ , VASP automatically extrapolates to the infinite basis set limit using a linear regression to the equation ^[3],^[2],^[4]:

$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 $\times$ ENCUTGW (again we do not recommend to change this default).

The integral over $\omega$ is evaluated by means of a highly accurate mini-max integration.^[5] The number of $\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 $\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.).

Large Systems

Virtually the same flags and procedures apply to the new low scaling RPA algorithm implemented in vasp.6.^[6]. However, in the last step ALGO=ACFDT needs to be replaced by either ALGO=ACFDTR or ALGO=RPAR in order to calculate the independent particle polarizability $\chi$ efficiently. To this end, the polarizability is calculated first on the imaginary time axis by contracting two non-interacting Green's functions Cite error: Closing </ref> missing for <ref> tag ^[2] ^[3] ^[4] ^[5] ^[6] ^[7] </references>

Contents

↑ Cite error: Invalid <ref> tag; no text was provided for refs named paier
↑ ^a ^b ^c J. Harl, L. Schimka, and G. Kresse, Phys. Rev. B 81, 115126 (2010).
↑ ^a ^b J. Harl and G. Kresse, Phys. Rev. B 77, 045136 (2008).
↑ ^a ^b J. Klimeš, M. Kaltak, and G. Kresse, Phys. Rev. B 90, 075125 (2014).
↑ ^a ^b M. Kaltak, J. Klimeš, and G. Kresse, Journal of Chemical Theory and Computation 10, 2498-2507 (2014).
↑ ^a ^b M. Kaltak, J. Klimeš, and G. Kresse, Phys. Rev. B 90, 054115 (2014)
↑ H. N. Rojas, R. W. Godby and R. J. Needs, Phys. Rev. Lett. 74, 1827 (1995)

[paier-1] Cite error: Invalid <ref> tag; no text was provided for refs named paier

[harl-2] J. Harl, L. Schimka, and G. Kresse, Phys. Rev. B 81, 115126 (2010).

[harl2-3] J. Harl and G. Kresse, Phys. Rev. B 77, 045136 (2008).

[klimes-4] J. Klimeš, M. Kaltak, and G. Kresse, Phys. Rev. B 90, 075125 (2014).

[kaltak-5] M. Kaltak, J. Klimeš, and G. Kresse, Journal of Chemical Theory and Computation 10, 2498-2507 (2014).

[kaltak2-6] M. Kaltak, J. Klimeš, and G. Kresse, Phys. Rev. B 90, 054115 (2014)

[rojas-7] H. N. Rojas, R. W. Godby and R. J. Needs, Phys. Rev. Lett. 74, 1827 (1995)

[1]

[2]

[3]

[4]

[5]

[6]

[7]