Kohn-Sham random-phase approximation

Electron correlation energies within the random-phase approximation can be calculated by the programs DIRPA, RPAX2 and ACFDT that are subdirectives of the driver command KSRPA. These methods should be used in conjunction with Kohn-Sham reference determinants, i.e., orbitals and orbital energies from a preceeding DFT calculation should be supplied.

A typical input to calculate the RPAX2 correlation energy is given by:

basis={
set,orbital; default,<basis>
set,mp2fit;  default,<basis>/mp2fit}

ks,pbe
{ksrpa; rpax2,orb=2100.2}

All methods are implemented using density-fitting of the two-electron repulsion integrals, see Refs. [2,4]. Because of this, auxiliary basis sets for fitting occupied-virtual orbital pairs have to be given, see section auxiliary basis sets for density fitting or resolution of the identity. No point-group symmetry can be used for all methods described in this section.

References:

RPA:
$[1]$ F. Furche, Phys. Rev. B 64, 195120 (2001).
RPAX2:
$[2]$ A. Heßelmann, Phys. Rev. A 85, 012517 (2012).
$[3]$ A. Heßelmann, Top. Curr. Chem. 365 97 (2015)
ACFDT(ALDA):
$[4]$ A. Heßelmann and A. Görling, J. Chem. Theory Comput. 9, 4382 (2013).

The direct RPA program (implemented with the algorithm described in [1]) has the following options:

  • ORB record number containing the orbital coefficients and eigenvalues (mandatory)
  • AUXBAS string containing the label for the auxiliary basis set (default MP2FIT)
  • CORE number of core orbitals (which are not correlated)
  • MAXIT maximum number of iterations (default ’40’)
  • THREN threshold for convergence of energy (default ’1d-8’)
  • FMIX mixing factor for the amplitude update $T^{\text{new}}=f T^{\text{old}}+(1-f)T^{\text{new}}$ (default: ’0.4d0’)
  • RESTART logical flag to enable a restart from an unfinished calculation. For this, the 3-index Coulomb integrals (Lc.dat) and the two amplitude files (T0.dat and T1.dat) are required if MODE=1 or MODE=2, see below. In case of MODE=3, Molpro’s file 4 needs to be saved in the previous calculation (using, e.g., file,4,rpa.dat at the beginning of the input file) (default: ’0’)
  • NOMAX maximum number of $N_{\text{aux}}\times N_{\text{virt}}$ batches to be kept in memory ($N_{\text{aux}}$: number of auxiliary basis functions, $N_{\text{virt}}$: number of virtual orbitals). (assumes: MODE=1, default: ’50’)
  • MODE can have the values ’0’,’1’,’2’ and ’3’. If MODE=0 all 3-index quantities are kept in memory. If MODE=1, external 3-index integral and amplitude files are written (to the wavefunction directory). With MODE=2 larger batches of amplitudes can be read/written as specified with NOMAX. MODE=3 is a duplicate of MODE=1, but the 3-index quantities are written to an internal Molpro file (default: ’MODE=3’)
  • L string containing the scratch file name for the 3-index Coulomb integrals (default: ’Lc.dat’)
  • T0 string containing the scratch file name for the amplitudes (default: ’T0.dat’)
  • T1 string containing the scratch file name for the amplitude updates (default: ’T1.dat’)
  • SOSEX logical flag, set to SOSEX=1 if the SOSEX energy shall be calculated after convergence (default: ’SOSEX=0’)

Note that in case of MODE=1 or MODE=2 it is recommended to have the wavefunction (wfu) directory located on a scratch partition. E.g., add the command line option -W /scratch/$USER/wfu to the Molpro command.

The RPAX2 method is an extension to the RPA and accounts for higher order particle-hole pair exchange contributions [2,3]. The RPAX2 program has the same options as the DIRPA program, see section DIRPA program. The following list shows a few additional ones that can be used:

  • DIR if set to DIR$\ne$0 this enables a direct RPA calculation (default ’0’)
  • MEM if set to MEM$\ne$0 the 3-index Coulomb integrals and the amplitudes are kept in memory (default: 0)

Spin-unrestricted calculations can be done using the URPAX2 program. In this case the orbitals from a preceeding unrestricted Kohn-Sham calculation have to be passed to the program (via the ORB key).

The ACFDT (adiabatic connection fluctuation-dissipation theorem) method is an alternative approach to derive the RPA. If used in conjunction with local adiabatic exchange-correlation kernels, the method can also describe electron-electron interaction contributions beyond the RPA. Currently, the ALDA xc-kernel can be used in the program (ACFDT(ALDA) method), see also Ref. [4]. The ACFDT program has the following options:

  • ORB record number containing the orbital coefficients and eigenvalues (mandatory)
  • AUXBAS string containing the label for the auxiliary basis set (default MP2FIT)
  • CORE number of core orbitals (which are not correlated)
  • NFREQ number of frequency quadrature points (default ’20’)
  • NCOUP number of coupling strength quadrature points (default ’7’)
  • GRIDTHR threshold for grid accuracy (default ’1d-10’)
  • THRKERN threshold for density in kernel integration (default ’1d-12’)
  • XFAC factor tor ALDA exchange contribution (default ’1d0’)
  • CFAC factor tor ALDA correlation contribution (default ’1d0’)
  • NOXC can be set to ’NOXC$\ne 0$’ if no ALDA xc-contribution shall be added to the electron-electron interaction (which then corresponds to a standard direct RPA calculation) (default ’0’)
  • OMQUAD set to ’1’ for Gauss-Chebyshev quadrature and ’2’ for Gauss-Legendre quadrature (default ’2’)
  • W0 parameter for Gauss-Legendre quadrature, see R. D. Amos et al., J. Phys. Chem. 89 (1985) 2186
  • L3ALPHA logical flag to switch on calculation of coupling-strength dependent xc-kernel integrals (deactivated by default)
  • FXC2IDX logical flag for performing a double density fitting approximation of the electron-electron interaction matrix (deactivated by default)
  • FXC2 a switch for various approaches to calculate the 2-index xc-kernel integrals (default ’1’)
  • THRDUM threshold for density on dummy centre quadrature points (if set to a large value, the dummy centre quadrature points are skipped completely) (default ’0d0’)
  • SING logical flag to enable handling of singularity of 2-index xc-kernel integrals, see also the following two options (enabled by default)
  • ESHIFT corresponds to $\epsilon$ in $f=1-\rho/(\rho+\epsilon)$ (default: ’1d-5’)
  • FSCAL corresponds to $s$ in $\rho=\rho+s\cdot f$ (default ’1d-4’)

The ACFDT(ALDA) method is ill-defined for short electron-electron distances, see Ref. [4] and F. Furche and T. Van Voorhis, J. Chem. Phys. 122, 164106 (2005). Because of this, the method does not have a defined basis set limit and the ACFDT program should not generally be used to calculate ACFDT(ALDA) correlation energies. Instead, the ACFDT2 should be utilised which implements the hybrid approach as described in Ref. [4] and which has, in addition to the ones given above, the following options:

  • SCAL scaling factor for the RPA kernel (used for short electron-electron distances) (default ’1d0’)
  • MU range-separation parameter. If not used, the program does not perform a correction for the short range electron-electron interaction.

For applying the correction as described in Ref. [4], the vaules of SCAL and MU have to be set to the values SCAL=0.6 and MU=2d0.

The ACFDT3 program implements an approximation to the ACFDT(ALDA) method assuming that the xc-kernel matrix depends linearly on the coupling strength (which is true for the exchange contribution but not, in general, for the correlation contribution to the kernel). Within this approximation the coupling-strength integration can be done analytically leading to a performance improvement over the ACFDT and ACFDT2 programs. The options for ACFDT3 are identical to the ones given above for ACFDT and ACFDT2.