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. An alterntive implementation of various RPA electron correlation methods by Julien Toulouse and co-workers is described in the subsection Random-phase approximation (RPATDDFT) program. The second alternative implementation of RPA is based on the adiabatic-connection fluctuation-dissipation and is described in the subsection RIRPA program. This implementation additionally allows calculations with σ-functionals, which at negligible computational cost significantly improves over RPA in the calculation of various chemical properties.
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:
MP2FIT)file,4,rpa.dat at the beginning of the input file) (default: ’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:
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:
MP2FIT)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:
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.
The random-phase approximation program (rpatddft) can be used to calculate RPA correlation energies after a SCF calculation. Additionnally, it can be used to calculate dynamic dipole polarizabilities, C$_6$ dispersion coefficients, and excitation energies. The program currently works without point-group symmetry.
List of the main keywords:
ECORR, <list of methods> Calculation of RPA correlation energies [1] (see options below)
PROPERTIES, <list of methods> Calculation of dynamic dipole polarizabilities and C$_6$ dispersion coefficients [10] (see options below)
EXCIT, <list of methods> Calculation of excitation energies [11] (see options below)
ORB,<orbrec> Record for input orbitals (required).as well as contextual options (see later for an explanation on this):
INTAC,NLAMBDA=<n>[,LAMBDA=<lambda>,WEIGHT=<weight>]
Number of quadrature points for the Gauss-Legendre numerical integration along the adiabatic connection for RPA calculations (default is 7). If LAMBDA and WEIGHT are given, assumes a one-point quadrature with given abscissa and weight.
INTFREQ,METHFREQ=<methfrq>,NFREQ=<nfreq>
Options for the numerical integration over the frequency variable of RPA calculations. METHFREQ governs the type of quadrature used (0(default) is Gauss-Chebyshev, 1 and 2 are Gauss-Legendre, 3 is Clenshaw-Curtis) and NFREQ governs the number of quadrature points (default is 16).
DIELMODE,mode Use the solid-state variant when performing DIEL calculations.
and global options, shared by all commands inside the rpatddft block:
STAB Check matrices stability conditions in RPA calculations. When used without an ECORR, EXCIT or PROPERTIES keyword, check the Hessian and RPA matrices eigenvalues and do nothing more.DFTKERNEL,<funcx>[,<funcc>]
Specify the exchange and correlation kernel for EXCIT (if only one argument is given, it is understood to be the exchange-correlation kernel).
C6 Computes C6 coefficients from last two saved polarizabilities.TDA Tamm-Dancoff approximation for EXCIT and PROPERTIES.NOMP2 The MP2 energy is calculated in certain situations where it is available almost for free, provided that some matrices are allocated. This behavior can be switched off by this NOMP2 keyword.NOSPINBLOCK For spin-unrestricted calculations, use a formalism where matrices are of $\alpha\alpha+\alpha\beta+\beta\alpha+\beta\beta$ dimensions (the default is to use a formalism with a nospinflip/splinflip block structure)NOSPINFLIP Exclude spin-flip dimensions of unrestricted RPA calculations that use the NOSPINBLOCK formalism (not suitable for all RPA variants).WRITEFILE Write files with eigenvalues, virtual orbital energies, dipole moments, dipole velocities, dipole accelerations and amplitudes from a TDA calculationVIAXPY #to comment#INTEGRAL,<nbr> Specify the two-electron integral transformation routine: 0 (still the default for spin-unrestricted and gradient calculations) is the ‘old’ one, 1 is the ‘old’ one that has been cleaned up, and 2 (default otherwise) is a much more efficient transformation using Molpro’s transform routine.OCC,<nocc> Explicitly specify the number of occupied orbitals (useful for fake pseudopotential calculations).CORE,<core> Specify core orbitals (default: last specified core orbitals or, if none, atomic inner shells)PRINT,<nbr> Level of print expected from the output (from 0(default) to 3).PRINT_INT,<nbr> Level of print of integrals (AO,MO,Orbtials,…), from 0 to 4.PRINT_TIME,<nbr> If greater or equal to 1, will print out information on time spent in routines.
Calculation of RPA correlation energies ECORR, <list of methods>
If no method is given, a SO2-RCCD calculation will be done (see below).
There are two main RPA variants [1]: dRPA (direct RPA, without the inclusion of an Hartree-Fock exchange kernel in the response function) and RPAx (with the Hartree-Fock exchange kernel included in the response function).
There are four main formulations in which the RPA equations can be derived. The adiabatic-connection fluctuation-dissipation theorem (ACFDT) equation involves integrations both over frequency and coupling constant: an analytical integration along the frequency variable followed by a quadrature on the coupling constant yields the adiabatic connection formulation (AC) [1], while an analytic integration on the coupling constant followed by a numerical integration over the frequency yields the dielectric formulation (DIEL) [2]. Two other formulations are obtained when the two integrations are carried out analytically: the plasmon formula (PLASMON) [1] and the ring coupled cluster doubles formulation (RCCD) [3]. These four formulations are not in general equivalent.
Most variants+formulations can readily by used in a spin-unrestricted context [6]. This is implemented in the code and does not need any further input from the user: the RPA program recognizes the spin-unrestricted character of a SCF calculation that was done beforehand and acts accordingly.
Gradients of most of the RCCD-formulation RPA energies are available, both without range-separation with RHF orbitals and with range-separation with RSH orbitals [9]. The calculations are triggered by the presence of the keyword FORCE or OPTG after the energy-related section (see examples at the end of the section).
The user can test the RPA program using make rpatddfttest, which proposes a variety of tests for RPA correlation energy calculations.
The keywords for the methods are constructed on the model:
<variant>-<formulation>-<alternative>
For the AC formulation, the available methods are:
DRPAI-AC dRPA calculation (see Refs. [1] and [4]).DRPAII-AC dRPA calculation, using antisymmetrized two-electron integrals (see Refs. [1] and [4]).RPAXI-AC RPAx calculation (see Refs. [1] and [5]).RPAXII-AC RPAx calculation, using antisymmetrized two-electron integrals (see Refs. [1] and [5]).DRPAI-AC-ALT dRPA calculation using an alternative derivation (see Ref. [7]).RPAXI-AC-ALT RPAx calculation using an alternative derivation (see Ref. [7]).DRPAI-AC-NOINT dRPA calculation without integration along the adiabatic connection (using the “kinetic” and “potential” contributions, sometimes called the “alternative plasmon formula”, see Ref. [1]).RPAXII-AC-NOINT RPAx calculation without integration along the adiabatic connection (using the “kinetic” and “potential” contributions, sometimes called the “alternative plasmon formula”, see Ref. [1]).For the DIEL formulation, the available methods are:
DRPAI-DIEL dRPA calculation (see Ref. [2]).DRPAIIA-DIEL dRPA-IIa approximation (see Ref. [2]).RPAXIA-DIEL RPAx-Ia approximation (see Ref. [2]).For the RCCD formulation, the available methods are:
DRPAI-RCCD dRPA-I calculation (see Ref. [3]).RPAXII-RCCD RPAx-II calculation (Szabo-Ostlund variant 1 is calculated too, see Ref. [3]).SO2-RCCD Szabo-Ostlund variant 2 (see Ref. [3]).SOSEX-RCCD dRPA-I+SOSEX correction (see Ref. [3]).RPAX2-RCCD RPAX2 approximation (see Ref. [8]).For the PLASMON formulation, the available methods are:
DRPAI-PLASMON dRPA-I calculation (see Ref. [1]).RPAXII-PLASMON RPAx-II calculation (see Ref. [1]).Note that to all these keywords are associated energy variables defined as :
ECORR_VARIANT_FORMULATION_ALTERNATIVE
(see the examples below).
Example of a dRPA-I calculation using the PBE functional:
{rks,pbe,orbital,2101.2}
{rhf,start=2101.2,maxit=0}
{rpatddft;
orb,2101.2;
ecorr,DRPAI-AC
}
e=ECORR_DRPAI_AC
Example of a range-separated RPAx-I calculation using the short-range PBE exchange-correlation functional and the range-separated parameter mu=0.5:
{int;erf,0.5}
{rks,exerfpbe,ecerfpbe;rangehybrid;orbital,2101.2}
{rpatddft;
orb,2101.2;
ecorr,RPAXI-AC
}
Example of several RPA calculations in the same run:
{rhf,orbital,2101.2}
{rpatddft;
orb,2101.2;
ecorr,DRPAI-AC,RPAXII-RCCD,DRPAI-DIEL
}
e1=ECORR_DRPAI_AC
e2=ECORR_RPAXII_RCCD
e3=ECORR_DRPAI_DIEL
(this way, the calculations are done with the same transformed integrals, i.e. without redoing the integral transformation).
Example of a dRPA-I gradient calculation:
{rks,pbe,orbital,2101.2}
{rpatddft;
orb,2101.2;
ecorr,DRPAI-RCCD
}
force
Example of a geometry optimization at the LDA+dRPA-I level:
{int;erf,0.5}
{rks,exerf,ecerf;rangehybrid;orbital,2101.2}
{rpatddft;
orb,2101.2;
ecorr,DRPAI-RCCD
}
optg
Calculation of properties, excitation energies and oscillator strengths
EXCIT, METHOD=<method> The EXCIT calculations output shows the excitation energies in ua, eV and nm, the oscillator strengths in length and velocity gauge, as well as the major excitations involved in each mode. The methods available are:
DRPA Direct random-phase approximation (or time-dependent Hartree).TDHF Time-dependent Hartree-Fock.TDDFT Time-dependent density-functional theory.RS-TDDFT Range-separated time-dependent density-functional theory [11].The exchange density functionals (FUNCX) available are:
LDAXERF (short-range LDA exchange density functional for the erf interaction [12]).The correlation density functionals (FUNCC) available are:
LDAC (Perdew-Wang-92 LDA correlation density functional)LDACERF (short-range LDA correlation density functional for the erf interaction [13]).Example of a range-separated time-dependent density-functional theory calculation using the short-range LDA exchange-correlation functional and the range-separated parameter mu=0.5:
mu=0.5
{int;erf,mu;save}
{rks,exerf,ecerf;rangehybrid;orbital,2101.2}
{int}
{setmu,mu}
{rpatddft;
orb,2101.2;
excit,method=rs-tddft;
dftkernel,funcx=ldaxerf,funcc=ldacerf
}
Example of a TDHF-TDA calculation with writing of several files for interfacing with a real-time propagation code (see Ref. [14]):
{rpatddft;
integral,1;
writefile,transmom.dat,energies.dat,virtual.dat,transmom_v.dat,amplitudes.dat,transmom_a.dat;
orb,2330.2;
excit,method=tdhf;
NOSPINBLOCK;
tda}
The files are: transmom.dat: transition moments in position form; energies.dat: total energies of the states; virtual.dat: virtual positive-energy orbitals; transmom_v.dat: transition moments in velocity form; amplitudes.dat: coefficients of excited-state wave functions over single-excited determinants; transmom_a.dat: transition moments in acceleration form
References
$[1]$ J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, J. Chem. Theory Comput. 7, 3116 (2011).
$[2]$ G. Jansen, B. Mussard, D. Rocca, J. G. Ángyán (in prep).
$[3]$ J. Toulouse, W. Zhu, A. Savin, G. Jansen, and J. G. Ángyán, J. Chem. Phys. 135, 084119 (2011).
$[4]$ F. Furche, Phys. Rev. B 64, 195120 (2001).
$[5]$ J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. Ángyán, Phys. Rev. Lett. 102, 096404 (2009).
$[6]$ B. Mussard, P. Reinhardt, J. G. Ángyán, and J. Toulouse, J. Chem. Phys. (submitted).
$[7]$ Heßelmann, A., Görling, A., Phys. Rev. Lett. 106, 093001 (2011).
$[8]$ Heßelmann, A., Phys. Rev. A 85, 012517 (2012).
$[9]$ B. Mussard, P. G. Szalay, J. G. Ángyán, J. Chem. Theory Comput. 10, 1968 (2014).
$[10]$ J. Toulouse, E. Rebolini, T. Gould, J. F. Dobson, P. Seal, J. G. Ángyán, J. Chem. Phys. 138, 194106 (2013).
$[11]$ E. Rebolini, A. Savin, J. Toulouse, Mol. Phys. 111, 1219 (2013).
$[12]$ J. Toulouse, A. Savin, and H.-J. Flad, Int. J. Quantum Chem. 100, 1047 (2004).
$[13]$ S. Paziani, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet, Phys. Rev. B 73, 155111 (2006).
$[14]$ E. Coccia, B. Mussard, M. Lebeye, J. Caillat, R. Taieb, J. Toulouse, and E. Luppi, Int. J. Quant. Chem. 116, 1120 (2016).
The RIRPA and URIRPA programs allow non-self-consistent spin-restricted and spin-unrestricted resolution of identity (RI) random phase approximation (RPA) [1-3] and σ-functional [4-6] calculations. These methods should be used in conjunction with conventional Kohn-Sham (KS) density functional theory (DFT) calculations, i.e. data from a preceding KS DFT calculation should be provided. Conventional KS DFT means calculations with LDA, GGA and hybrid exchange-correlation functionals.
Bibilography:
RPA:
[1] F. Furche, Phys. Rev. B, 195120 (2001)
[2] A. Heßelmann and A. Görling, Mol. Phys. 109, 2473 (2011)
[3] X. Ren, P. Rinke, C. Joas, and M. Scheffler, J. Mater. Sci. 47, 7447 (2012)
σ-functionals:
[4] E. Trushin, A. Thierbach, A. Görling, J. Chem. Phys. 154, 014104 (2021)
[5] S. Fauser, E. Trushin, C. Neiss, A. Görling, J. Chem. Phys. 155, 134111 (2021)
[6] C. Neiss, S. Fauser, A. Görling, J. Chem. Phys. 158, 044107 (2023)
Other papers cited in the documentation:
[7] E. Trushin, A. Görling, J. Chem. Phys. 155, 054109 (2021)
We kindly ask you to cite the original publications of the corresponding methods in the publications that result from these programs.
Example input file for spin-restricted calculations for the CO molecule:
gthresh,twoint=1d-20,energy=1d-10,orbital=1d-8 ! tighter thresholds are recommended
gdirect ! integral-direct mode
basis={
default,aug-cc-pwCVQZ ! orbital basis
set,ri; default,aug-cc-pwCVQZ/mp2fit ! RI basis
}
symmetry,nosym ! RIRPA does not use symmetry
angstrom
geometry={
2
C 0.000000 0.000000 -0.646514
O 0.000000 0.000000 0.484886
}
df-ks,pbe,df_basis=aug-cc-pwCV5Z/mp2fit,maxit=200 ! DFT calculation
{cfit,basis_coul=aug-cc-pwCV5Z/mp2fit,basis_exch=aug-cc-pwCV5Z/mp2fit}
acfd;rirpa ! RPA/sigma-functional calculation; one can use alternatively: ksrpa;rirpa
As well as an example of spin-unrestricted calculation for the NH2 molecule:
gthresh,twoint=1d-20,energy=1d-10,orbital=1d-8 ! tighter thresholds are recommended
gdirect ! integral-direct mode
basis={
default,aug-cc-pwCVQZ ! orbital basis
set,ri; default,aug-cc-pwCVQZ/mp2fit ! RI basis
}
symmetry,nosym ! URIRPA does not use symmetry
angstrom
geometry={
3
N 0.000000 0.000000 0.142235
H 0.000000 0.800646 -0.497821
H 0.000000 -0.800646 -0.497821
}
spin=1
df-uks,pbe,df_basis=aug-cc-pwCV5Z/mp2fit,maxit=200 ! DFT calculation
{cfit,basis_coul=aug-cc-pwCV5Z/mp2fit,basis_exch=aug-cc-pwCV5Z/mp2fit}
acfd;urirpa ! RPA/sigma-functional calculation; one can use alternatively: ksrpa;urirpa
The following options are available for the RIRPA and URIRPA programs: