Correlation energies from TDDFT

The Random Phase Approximation (RPA) for correlation energies comprises a nice non-empirical expression for the correlation energy that can be naturally combined with exact exchange to calculate binding energies. Due to the non-local nature of the approximation, RPA gives a good account of dispersive forces and is the only xc approximation capable of describing intricate binding regimes where covalent and van der Waals interactions are equally important. However, RPA does not provide a very accurate description of strong covalent bonds and typically performs worse than standard GGA functionals for molecular atomization energies and cohesive energies of solids.

In general, the exact correlation energy can be expressed in terms of the exact response function as

\[E_c = -\int_0^{\infty}\frac{d\omega}{2\pi}\int_0^1d\lambda\text{Tr}\Big[v\chi^\lambda(\omega)-v\chi^{KS}(\omega)\Big],\]

and the RPA approximation for correlation energies is simply obtained from the RPA approximation for the response function. From the point of view of TDDFT, the response function can be expressed exactly in terms of the Kohn-Sham response function and the exchange-correlation kernel \(f_{xc}\):

\[\chi^\lambda(\omega) = \chi^{KS}(\omega) + \chi^{KS}(\omega)\Big[\lambda v+f_{xc}^\lambda(\omega)\Big]\chi^\lambda(\omega).\]

The RPA is obtained by neglecting the exchange-correlation kernel and it should be possible to improve RPA by including simple approximations for this kernel. However, if one tries to use a simple adiabatic kernel, one encounters severe convergence problems and results become worse than RPA! The reason for this is that the locality of adiabatic kernels renders the pair-correlation function divergent. As it turns out, the adiabatic correlation hole can be renormalized by a simple non-empirical procedure, which results in a density-dependent non-locality in the kernel. This can be done for any adiabatic kernel and the method has implemented for LDA and PBE. We refer to these approximations as renormalized adiabatic LDA (rALDA) and renormalized adiabatic PBE (rAPBE). We only include the exchange part of the kernel, since this part is linear in \(\lambda\) and the kernel thus only needs to be evaluated for \(\lambda=1\).

For more details on the theory and implementation of RPA we refer to RPA correlation energy and the tutorial Calculating RPA correlation energies. The RPA tutorial should be studied before the present tutorial, which inherits much of the terminology from RPA. Details on the theory, implementation and benchmarking of the renormalized kernels can be found in Refs. [1], [2], and [3].

Below we give examples on how to calculate the correlation energy of a Hydrogen atom as well as the rAPBE atomization energy of a \(CO\) molecule and the rAPBE cohesive energy of diamond. Note that some of the calculations in this tutorial will need a lot of CPU time and are essentially not possible without a supercomputer.

Finally, we note that there is some freedom in deciding how to include the density dependence in the kernel. By default the kernel is constructed from a two-point average of the density. However as shown in Example 4 it is possible to instead use a reciprocal-space averaging procedure. Within this averaging scheme it is possible to explore different approximations for \(f_{xc}\), for instance a simple dynamical kernel, or a jellium-with-gap model, which displays a \(1/q^2\) divergence for small \(q\). More details can be found in [4].

Example 1: Correlation energy of the Hydrogen atom

As a first demonstration of the deficiencies of RPA, we calculate the correlation energy of a Hydrogen atom. The exact correlation energy should vanish for any one-electron system. The calculations can be performed with the following scripts, starting with a standard DFT-LDA calculation:

from ase import Atoms
from ase.parallel import paropen
from gpaw import GPAW
from gpaw import PW

resultfile = paropen('H.ralda.DFT_corr_energies.txt', 'w')
resultfile.write('DFT Correlation energies for H atom\n')

H = Atoms('H', [(0, 0, 0)])
H.center(vacuum=2.0)
calc = GPAW(mode=PW(400, force_complex_dtype=True),
            parallel={'domain': 1},
            hund=True,
            txt='H.ralda_01_lda.output.txt',
            xc='LDA')

H.calc = calc
E_lda = H.get_potential_energy()
E_c_lda = -calc.get_xc_difference('LDA_X')

resultfile.write(f'LDA correlation: {E_c_lda} eV')
resultfile.write('\n')

calc.diagonalize_full_hamiltonian()
calc.write('H.ralda.lda_wfcs.gpw', mode='all')

followed by an RPA calculation:

from gpaw.xc.rpa import RPACorrelation

rpa = RPACorrelation('H.ralda.lda_wfcs.gpw',
                     ecut=300,
                     txt='H.ralda_02_rpa_at_lda.output.txt')
rpa.calculate()

and finally one using the rALDA kernel:

from gpaw.xc.fxc import FXCCorrelation

fxc = FXCCorrelation('H.ralda.lda_wfcs.gpw',
                     xc='rALDA', txt='H.ralda_03_ralda.output.txt',
                     ecut=300)
fxc.calculate()

The analogous set of scripts for PBE/rAPBE are H.ralda_04_pbe.py, H.ralda_05_rpa_at_pbe.py and H.ralda_06_rapbe.py. The computationally heavy RPA/rALDA/rAPBE parts can be parallelized efficiently on multiple CPUs. After running the scripts the LDA and PBE correlation energies may be found in the file H.ralda.DFT_corr_energies.txt. Note that a rather small unit cell is used and the results may not be completely converged with respect to cutoff and unit cell. Also note that the correlation energy is calculated at different cutoff energies up to 300 eV and the values based on two-point extrapolation are printed at the end (see Calculating RPA correlation energies and RPA correlation energy for a discussion on extrapolation). The results in eV are summarized below.

LDA

RPA

rALDA

-0.56

-0.56

-0.032

PBE

RPA

rAPBE

-0.15

-0.56

-0.009

The fact that RPA gives such a dramatic underestimation of the correlation energy is a general problem with the method, which is also seen for bulk systems. For example, for the homogeneous electron gas RPA underestimates the correlation energy by ~0.5 eV per electron for a wide range of densities.

Example 2: Atomization energy of CO

Although RPA severely underestimates absolute correlation energies in general, energy differences are often of decent quality due to extended error cancellation. Nevertheless, RPA tends to underbind and performs slightly worse than PBE for atomization energies of molecules. The following example shows that rAPBE not only corrects the absolute correlation energies, but also performs better than RPA for atomization energies.

First we set up a ground state calculation with lots of unoccupied bands. This is done with the script:

from ase import Atoms
from ase.parallel import paropen
from gpaw import GPAW, FermiDirac
from gpaw.mixer import MixerSum
from gpaw.hybrids.energy import non_self_consistent_energy as nsc_energy
from gpaw import PW

# CO

CO = Atoms('CO', [(0, 0, 0), (0, 0, 1.1283)])
CO.center(vacuum=2.7)
calc = GPAW(mode=PW(600, force_complex_dtype=True),
            symmetry='off',
            parallel={'domain': 1},
            xc='PBE',
            txt='CO.ralda_01_CO_pbe.txt',
            convergence={'density': 1.e-6})

CO.calc = calc
E0_pbe = CO.get_potential_energy()

E0_hf = nsc_energy(calc, 'EXX').sum()

calc.diagonalize_full_hamiltonian()
calc.write('CO.ralda.pbe_wfcs_CO.gpw', mode='all')

# C

C = Atoms('C')
C.set_cell(CO.cell)
C.center()
calc = GPAW(mode=PW(600, force_complex_dtype=True),
            symmetry='off',
            parallel={'domain': 1},
            xc='PBE',
            mixer=MixerSum(beta=0.1, nmaxold=5, weight=50.0),
            hund=True,
            occupations=FermiDirac(0.01),  # , fixmagmom=True),
            txt='CO.ralda_01_C_pbe.txt',
            convergence={'density': 1.e-6})

C.calc = calc
E1_pbe = C.get_potential_energy()

E1_hf = nsc_energy(calc, 'EXX').sum()

f = paropen('CO.ralda.PBE_HF_C.dat', 'w')
print(E1_pbe, E1_hf, file=f)
f.close()

calc.diagonalize_full_hamiltonian()
calc.write('CO.ralda.pbe_wfcs_C.gpw', mode='all')

# O

O = Atoms('O')
O.set_cell(CO.cell)
O.center()
calc = GPAW(mode=PW(600, force_complex_dtype=True),
            symmetry='off',
            parallel={'domain': 1},
            xc='PBE',
            mixer=MixerSum(beta=0.1, nmaxold=5, weight=50.0),
            hund=True,
            txt='CO.ralda_01_O_pbe.txt',
            convergence={'density': 1.e-6})

O.calc = calc
E2_pbe = O.get_potential_energy()

E2_hf = nsc_energy(calc, 'EXX').sum()

calc.diagonalize_full_hamiltonian()
calc.write('CO.ralda.pbe_wfcs_O.gpw', mode='all')

f = paropen('CO.ralda.PBE_HF_CO.dat', 'w')
print('PBE: ', E0_pbe - E1_pbe - E2_pbe, file=f)
print('HF: ', E0_hf - E1_hf - E2_hf, file=f)
f.close()

which takes on the order of 6-7 CPU hours. The script generates three gpw files containing the wavefunctions, which are the input to the rAPBE calculation. The PBE and non self-consistent Hartree-Fock atomization energies are also calculated and written to the file CO.ralda.PBE_HF_CO.dat. Be aware that using symmetries (i. e. not using symmetry='off' in the calculator) may cause problems if you want to calculate non self-consistent HF energies for atoms and molecules. Next we calculate the RPA and rAPBE energies for CO with the script

from ase.parallel import paropen
from ase.units import Hartree
from gpaw.xc.rpa import RPACorrelation
from gpaw.xc.fxc import FXCCorrelation

fxc0 = FXCCorrelation('CO.ralda.pbe_wfcs_CO.gpw',
                      xc='rAPBE',
                      ecut=400,
                      nblocks=8,
                      txt='CO.ralda_02_CO_rapbe.txt')

E0_i = fxc0.calculate()

f = paropen('CO.ralda_rapbe_CO.dat', 'w')
for ecut, E0 in zip(fxc0.rpa.ecut_i, E0_i):
    print(ecut * Hartree, E0, file=f)
f.close()

rpa0 = RPACorrelation('CO.ralda.pbe_wfcs_CO.gpw',
                      ecut=400,
                      nblocks=8,
                      txt='CO.ralda_02_CO_rpa.txt')

E0_i = rpa0.calculate()

f = paropen('CO.ralda_rpa_CO.dat', 'w')
for ecut, E0 in zip(rpa0.ecut_i, E0_i):
    print(ecut * Hartree, E0, file=f)
f.close()

The energies for C and O are obtained from the corresponding scripts CO.ralda_03_C_rapbe.py and CO.ralda_04_O_rapbe.py. The results for various cutoffs are written to the files like CO.ralda_rpa_CO.dat and CO.ralda_rapbe_CO.dat. We also print the correlation energies of the C atom to be used in a tutorial below. As in the case of RPA the converged result is obtained by extrapolation using the script

import numpy as np
from gpaw.utilities.extrapolate import extrapolate

E_pbe, E_hf = np.genfromtxt('CO.ralda.PBE_HF_CO.dat')[:, 1]
assert abs(E_pbe - -11.74) < 0.01
assert abs(E_hf - -7.37) < 0.01

CO_rpa = np.loadtxt('CO.ralda_rpa_CO.dat')
C_rpa = np.loadtxt('CO.ralda_rpa_C.dat')
O_rpa = np.loadtxt('CO.ralda_rpa_O.dat')

a = CO_rpa
a[:, 1] -= (C_rpa[:, 1] + O_rpa[:, 1])
ext, A, B, sigma = extrapolate(a[:, 0], a[:, 1], reg=3, plot=False)
assert abs(A - -3.27) < 0.02, A
assert abs(E_hf + A - -10.64) < 0.01

if 0:
    # RAPBE not currently working!
    CO_rapbe = np.loadtxt('CO.ralda_rapbe_CO.dat')
    C_rapbe = np.loadtxt('CO.ralda_rapbe_C.dat')
    O_rapbe = np.loadtxt('CO.ralda_rapbe_O.dat')

    a = CO_rapbe
    a[:, 1] -= (C_rapbe[:, 1] + O_rapbe[:, 1])
    ext, A, B, sigma = extrapolate(a[:, 0], a[:, 1], reg=3, plot=False)
    assert abs(A - -3.51) < 0.01
    assert abs(E_hf + A - -10.88) < 0.01

If pylab is installed, the plot=False can be change to plot=True to visualize the quality of the extrapolation. The final results are displayed below

PBE

HF

RPA

rAPBE

Experimental

11.74

7.37

10.64

10.88

11.23

Example 3: Cohesive energy of diamond

The error cancellation in RPA works best when comparing systems with similar electronic structure. In the case of cohesive energies of solids where the bulk energy is compared to the energy of isolated atoms, RPA becomes even worse than for atomization energies of molecules. Here we illustrate this for the cohesive energy of diamond and show that the rAPBE approximation corrects the problem. The initial orbitals are obtained with the script

from ase import Atoms
from ase.build import bulk
from ase.dft import monkhorst_pack
from ase.parallel import paropen
from gpaw import GPAW, FermiDirac
from gpaw import PW
from gpaw.hybrids.energy import non_self_consistent_energy as nsc_energy
import numpy as np

# Monkhorst-Pack grid shifted to be gamma centered
k = 8
kpts = monkhorst_pack([k, k, k])
kpts += [1. / (2 * k), 1. / (2 * k), 1. / (2 * k)]

cell = bulk('C', 'fcc', a=3.553).get_cell()
a = Atoms('C2', cell=cell, pbc=True,
          scaled_positions=((0, 0, 0), (0.25, 0.25, 0.25)))

calc = GPAW(mode=PW(600),
            xc='PBE',
            occupations=FermiDirac(width=0.01),
            convergence={'density': 1.e-6},
            kpts=kpts,
            parallel={'domain': 1},
            txt='diamond.ralda_01_pbe.txt')

a.calc = calc
E_pbe = a.get_potential_energy()
E_hf = nsc_energy(calc, 'EXX').sum()
E_C = np.loadtxt('CO.ralda.PBE_HF_C.dat')

f = paropen('diamond.ralda.PBE_HF_diamond.dat', 'w')
print('PBE: ', E_pbe / 2 - E_C[0], file=f)
print('HF: ', E_hf / 2 - E_C[1], file=f)
f.close()

calc.diagonalize_full_hamiltonian()
calc.write('diamond.ralda.pbe_wfcs.gpw', mode='all')

which takes roughly 5 minutes. The script generates diamond.ralda.pbe_wfcs.gpw and uses a previous calculation of the C atom to calculate the EXX and PBE cohesive energies that are written to diamond.ralda.PBE_HF_diamond.dat. The RPA and rAPBE correlation energies are obtained with the script:

from ase.parallel import paropen
from ase.units import Hartree
from gpaw.xc.rpa import RPACorrelation
from gpaw.xc.fxc import FXCCorrelation

fxc = FXCCorrelation('diamond.ralda.pbe_wfcs.gpw', xc='rAPBE',
                     ecut=400,
                     nblocks=8,
                     txt='diamond.ralda_02_rapbe.txt')
E_i = fxc.calculate()

f = paropen('diamond.ralda.rapbe.dat', 'w')
for ecut, E in zip(fxc.rpa.ecut_i, E_i):
    print(ecut * Hartree, E, file=f)
f.close()

rpa = RPACorrelation('diamond.ralda.pbe_wfcs.gpw',
                     ecut=400,
                     nblocks=8,
                     txt='diamond.ralda_02_rpa.txt')
E_i = rpa.calculate()

f = paropen('diamond.ralda.rpa.dat', 'w')
for ecut, E in zip(rpa.ecut_i, E_i):
    print(ecut * Hartree, E, file=f)
f.close()

This takes on the order of 30 CPU hours, but can be parallelized efficiently. Finally the correlation part of the cohesive energies are obtained by extrapolation with the script

from gpaw.utilities.extrapolate import extrapolate
import numpy as np

E_pbe, E_hf = np.genfromtxt('diamond.ralda.PBE_HF_diamond.dat')[:, 1]
assert abs(E_pbe - -7.75) < 0.01
assert abs(E_hf - -5.17) < 0.01

a = np.loadtxt('diamond.ralda.rpa.dat')
b = np.loadtxt('CO.ralda_rpa_C.dat')
ext, A, B, sigma = extrapolate(a[:, 0], a[:, 1] / 2 - b[:, 1],
                               reg=3, plot=False)
assert abs(A - -1.89) < 0.01, A
assert abs(E_hf + A - -7.06) < 0.01

if 0:
    # RAPBE currently not working!
    a = np.loadtxt('diamond.ralda.rapbe.dat')
    b = np.loadtxt('CO.ralda_rapbe_C.dat')
    ext, A, B, sigma = extrapolate(a[:, 0], a[:, 1] / 2 - b[:, 1],
                                   reg=3, plot=False)
    assert abs(A - -1.48) < 0.01
    assert abs(E_hf + A - -6.65) < 0.01

The results are summarized below

PBE

HF

RPA

rAPBE

Experimental

7.75

5.17

7.06

6.65

7.55

As anticipated, RPA severely underestimates the cohesive energy, while PBE performs much better, and rAPBE comes very close to the experimental value.

Example 4: Correlation energy of diamond with different kernels

Finally we provide an example where we try different exchange-correlation kernels. For illustrative purposes we use a basic computational setup - as a result these numbers should not be considered converged! As usual we start with a ground state calculation to get the electron density and wavefunctions for diamond:

from ase.build import bulk
from gpaw import GPAW, FermiDirac
from gpaw import PW

bulk_c = bulk('C', a=3.5454859)
calc = GPAW(mode=PW(600.0),
            parallel={'domain': 1},
            xc='LDA',
            occupations=FermiDirac(width=0.01),
            kpts={'size': (6, 6, 6), 'gamma': True},
            txt='diam_kern.ralda_01_lda.txt')

bulk_c.calc = calc
E_lda = bulk_c.get_potential_energy()
calc.diagonalize_full_hamiltonian()
calc.write('diam_kern.ralda.lda_wfcs.gpw', mode='all')

The default method of constructing the kernel is to use a two-point density average. Therefore the following simple script gets the rALDA correlation energy within this averaging scheme:

from gpaw.xc.fxc import FXCCorrelation
from ase.parallel import paropen

fxc = FXCCorrelation('diam_kern.ralda.lda_wfcs.gpw', xc='rALDA',
                     ecut=[131.072],
                     txt='diam_kern.ralda_02_ralda_dens.txt',
                     avg_scheme='density')
E_i = fxc.calculate()

resultfile = paropen('diam_kern.ralda_kernel_comparison.dat', 'w')
resultfile.write(str(E_i[-1]) + '\n')
resultfile.close()

However, an alternative method of constructing the kernel is to work in reciprocal space, and average over wavevectors rather than density. To use this averaging scheme, we add the flag avg_scheme='wavevector':

from gpaw.xc.fxc import FXCCorrelation
from ase.parallel import paropen

fxc = FXCCorrelation('diam_kern.ralda.lda_wfcs.gpw', xc='rALDA',
                     ecut=[131.072],
                     txt='diam_kern.ralda_03_ralda_wave.txt',
                     avg_scheme='wavevector')
E_i = fxc.calculate()

resultfile = paropen('diam_kern.ralda_kernel_comparison.dat', 'a')
resultfile.write(str(E_i[-1]) + '\n')
resultfile.close()

Using this averaging scheme opens a few more possible choices for the kernel.

For further comparison we include an RPA calculation as well.

from gpaw.xc.fxc import FXCCorrelation
from ase.parallel import paropen

fxc = FXCCorrelation('diam_kern.ralda.lda_wfcs.gpw', xc='RPA',
                     ecut=[131.072, 80.0],
                     txt='diam_kern.ralda_08_rpa.txt')
E_i = fxc.calculate()

resultfile = paropen('diam_kern.ralda_kernel_comparison.dat', 'a')
resultfile.write(str(E_i[-1]) + '\n')
resultfile.close()

Here we summarize the above calculations and show the correlation energy/electron (in eV), obtained at an (unconverged) cutoff of 131 eV:

rALDA (dens. av.)

rALDA (wave. av)

RPA

-1.16

-1.13

-1.40

Incidentally, a fully converged RPA calculation gives a correlation energy of -1.781 eV per electron.

We conclude with some practical points. The wavevector-averaging scheme is less intuitive than the density-average, but avoids some difficulties such as having to describe the \(1/r\) divergence of the Coulomb interaction in real space. It also can be faster to construct for systems with many k points. Finally we point out that the rALDA and rAPBE kernels are also special because they have explicit spin-polarized forms.