RPA correlation energy¶
The correlation energy within the Random Phase Approximation (RPA) can be written
where \(\chi^0(i\omega)\) is the noninteracting (KohnSham) response function evaluated at complex frequencies, \(\text{Tr}\) is the Trace and \(\it{v}\) is the Coulomb interaction. The response function and Coulomb interaction are evaluated in a plane wave basis as described in Linear dielectric response of an extended system and Linear dielectric response of an extended system: theory and for periodic systems the Trace therefore involves a summation over \(\mathbf{q}\)points, which are determined from the Brillouin zone sampling used when calculating \(\chi^0(i\omega)\).
The RPA correlation energy is obtained by:
from gpaw.xc.rpa import RPACorrelation
rpa = RPACorrelation(calc, txt='rpa_correlation.txt')
E_rpa = rpa.calculate(ecut=400)
where calc is either a calculator object containing converged wavefunctions from a ground state calculation or a string reference to a .gpw file containing wavefunctions. If calc is a calculator object it should be loaded in serial since the RPA parallellization scheme is rather different from that of standard DFT calculatons. txt denotes the output file. The RPACorrelation also takes a number of optional keywords described below. The calculate() function performs the actual calculation at the cutoff energy specified by ecut (in eV). In addition the rpa calculator will calculate the correlation energy at four values for the cutoff energies up to the specified cutoff, but one can also give a list of cutoff values instead. By default, the response function is calculated with the same number of bands as the number of plane waves, but one can also specify that it should use N bands with nbands=N in the calculate() function.
Parameters¶
keyword  type  default value  description 

nfrequencies 
int 
16  Number of Gausslegendre points used in the integration. 
frequency_cut 
float 

The maximum frequency is the largest frequency included in the GaussLegendre integration. The integral is always an approximation to the infinite integral, but the max frequency determines the distribution of frequencies. 
frequency_scale 
float 
2.0 (eV)  The frequency scale sets the density of frequency points near \(\omega = 0\). 
frequencies 
numpy.ndarray 
None  Specifies frequency points used to integrate the correlation integrand. Ex: numpy.linspace(0,20,201). If None, the Gausslegendre method is used. 
weights 
numpy.ndarray 
None  Should be used in conjunction with frequencies (e.i. when not using the GaussLegendre integration). For example np.array([0.5,1,1,…,1,1,0.5] gives a trapezoid integration 
skip_gamma 
bool 
False  For metals the \(\mathbf{q} = 0\) point can give rise to divergent contributions and it may be faster to converge the kpoint sampling if this point is excluded. 
nblocks 
int 
1  Gvector parallelization. Default parallelization scheme is over kpoints, spin and bands. If memory becomes an issue it can be an advantage to use Gvector parallelization also. 
filename 
str 
None  Restart file. If calculations with kpoint sampling, the contributions from different qpoints are calculated sequentially and written to filename such that these do not have to be recalculated when a calculation is restarted. 
In addition to the usual kpoint and plane wave cutoff, the RPA correlation energy needs to be converged with respect to a plane wave cutoff in the response function (set by ecut) and the frequency integration. As it turns out, the integrand is usually rather smooth and one can perform the integration with 816 (special!) GaussLegendre frequency points, but see the tutorial Calculating RPA correlation energies for an example of converging the frequency integration.
Convergence¶
A major complication with the RPA correlation energy is that it converges very slowly with the number of unoccupied bands included in the evaluation of \(\chi^0(i\omega)\). However, as described in Ref. [1] the high energy part of the response function resembles the Lindhard function, which for high energies gives a correlation energy converging as
where \(E^{\chi}_{cut}\) is cutoff energy used in the evaluation of \(\chi^0\). With an external potential, the number of unoccupied bands is an additional convergence parameter, but for reproducing the scaling of the Lindhard function, it is natural to set the total number of bands equal to the number of plane waves used. Thus, to obtain a converged RPA correlation energy one should proceed in three steps.
 Perform a ground state calculation with a lot of converged unoccupied bands.
 Define a list of cutoff energies  typically something like [200, 225, 250, 275, 300] (eV). For each cutoff energy perform an RPA correlation energy calculation with the number bands \(n\) set equal to the number of plane waves defined by that cutoff energy.
 Fit the list of obtained correlation energies to \(E_c^{RPA}(E) = E_c^{\infty}+A/E^{3/2}\) to obtain \(E_c^{\infty}=E_c^{RPA}\).
Per default, the rpa module defines a list of five cutoff energies up to the specified value and performs the extrapolation at the end of the calculation. If one is not interested in the total correlation energy, but only energy differences between similar systems, it is sometimes possible to avoid the extrapolation procedure and the rpa correlation energy can be obtained at a single point by specifying a list with one element (for example ecut=[400]).
[1]  J. Harl and G. Kresse, Phys. Rev. B 77, 045136 (2008) 
[2]  J. Harl and L. Schimka and G. Kresse, Phys. Rev. B 81, 115126 (2010) 