RPA, GW and excited states

Methods built on the density–density response function: RPA correlation energies, GW quasiparticle energies, BSE and TDDFT excitations, and polarizabilities and \( C_6 \) coefficients. All need an RI auxiliary basis ([rpa] auxbasis). The [rpa], [gw] and [tddft] keys are in Input file; grades are on Capabilities and validation.

What PDEP does in ferric

The independent-particle response \( \chi_0(i\omega) \) is built in the RI auxiliary basis from the three-centre B tensors, as an explicit (Adler–Wiser) sum over every occupied–virtual pair \( ia \), weighted by \( 4\varepsilon_{ia}/(\omega^2 + \varepsilon_{ia}^2) \) (sternheimer::dielectric_matrix). The sum over empty states is still there.

PDEP (projective dielectric eigenpotentials) then works in the eigenbasis of the static dielectric matrix in that RI space. Eigenpotentials whose eigenvalue is within trunc_thresh of 1 (default 1e-4) carry almost no screening and are dropped, so the frequency-dependent work runs in a smaller basis. That compression, not the removal of the empty-state sum, is what PDEP contributes here. How much it saves depends on the threshold; runs that need the full-rank answer set trunc_thresh = 0.0, as the GW, BSE and \( C_6 \) examples do.

RPA correlation energy

What it is. Direct RPA (dRPA) correlation from the dielectric eigenvalues on an imaginary-frequency quadrature.

• PDEP-RPA, closed shell: method.kind = "pdep-rpa" (examples/water-pdep-rpa.toml); Python ferric.run_pdep_rpa. Proven.
• U-PDEP-RPA, open shell over a spin-summed dielectric. From the CLI, set method.kind = "pdep-rpa" with multiplicity > 1 and task = "energy": the CLI solves UHF (UKS with [rpa] xc) with MOM after 5 iterations and runs U-PDEP-RPA on it. It is CLI-only: Python run_pdep_rpa is closed shell only. The library (ferric_rpa::run_u_pdep_rpa) also accepts a ROHF (or ROKS) reference, which it semi-canonicalizes first: each spin uses the orbitals and energies of its own Fock matrix, diagonalized in its occupied and virtual blocks (see the anchors).
• Attenuated RPA: short-range correlation with an erfc operator.
• RS-MP2 + LR-RPA: short-range MP2 plus long-range dRPA, on the MP2 page.

The static eigensolve defaults to Lanczos, with a dense path for small problems. Geometry optimization with pdep-rpa is supported (task = "optimize") on a closed-shell RHF reference.

Exact and local dRPA

What it is. method.kind = "drpa" (Python ferric.run_drpa) is dRPA@HF computed by the drCCD Riccati equations in the Boys-localized basis, closed shell, energy only. It is exact by default: no amplitude is truncated (examples/water-drpa.toml), and the energy equals the canonical plasmon formula to ≤ 1e-12 Ha. Riccati, plasmon and full-rank PDEP-RPA are algorithms for the same exact dRPA; measured agreement is ≤ 2.6e-14 Ha (H2/STO-3G, water/6-31G with and without frozen core, PDEP with trunc_thresh = 0 at 64 Gauss–Legendre points). Proven (narrow).

Which exact algorithm. The Riccati solve holds a ring-product plan of no³·nv² numbers, no times the size of the amplitudes, so it is the small-system path: C12 thrashed and was then killed for memory. A run that cannot fit the memory budget is refused before the SCF and pointed at pdep-rpa with [rpa] trunc_thresh = 0, which gives the same energy (to its frequency-quadrature error) at far lower memory and is faster at every size measured (n-alkanes C4–C16).

The local approximation ([local] scheme = "amplitude-threshold" with eps, examples/water-drpa-local.toml; Python run_drpa(..., local="amplitude-threshold", eps=1e-4)) drops pair amplitudes whose localized |2(ia|jb)| is at or below eps. eps has no default and is printed and logged with the kept fraction; eps = 0 is the exact method. dRPA is not variational, so the error is first order in what is dropped. It is positive (less correlation) at every point measured, which is a measurement, not a guarantee, and grows faster than linearly in eps: on n-octane / 6-31G it is 6.7e-7, 3.9e-5, 5.5e-4 and 1.0e-2 Ha at eps = 1e-6, 1e-5, 1e-4 and 1e-3, keeping 82%, 50%, 18% and 3% of the amplitudes. [local] reference = true (Python compute_reference=True) also computes the canonical plasmon dRPA and prints the error against it. [local] eps_sweep (Python run_drpa_scan) evaluates several eps on one SCF and one localized assembly. Proven (narrow) through its exact limit.

What is measured (n-alkanes, 6-31G / cc-pVDZ-RI, frozen carbon cores, Coulomb, single thread; local error against the plasmon formula, PDEP error against full-rank PDEP at 32 quadrature points). The fraction of the truncated object retained at 1 kcal/mol error:

The amplitude threshold compresses more as the molecule grows, and PDEP's fraction stays flat. The two fractions are of different objects (no²nv² amplitudes against naux modes), so this is not a cost comparison. In wall time at a matched error of about 1 kcal/mol (reference off), truncated PDEP stays about 8× faster from C4 to C16 (C16: 165.7 s local at eps = 1e-4 against about 21 s PDEP), and both grow at the same rate in that range. No speedup over PDEP is claimed for local dRPA.

GW

What it is. Quasiparticle energies from the GW self-energy: G0W0, COHSEX, evGW0 and evGW, closed shell, plus unrestricted U-GW. The starting point is HF by default or a KS functional ([rpa] xc). From a KS starting point the static term Σx − v_xc enters the G0W0, evGW₀ and evGW quasiparticle equation (for U-GW, each spin's own equation with that spin's v_xc), so Σc is evaluated at the shifted root; COHSEX and U-COHSEX, which are static, add it to the quasiparticle energy.

Run it. method.kind = "gw" with [gw] method = "g0w0" (examples/water-g0w0-pbe.toml, open shell examples/oh-ugw.toml); Python ferric.run_gw, run_u_gw. The open-shell reference is UHF by default; [gw] reference = "rohf" (Python run_u_gw(reference="rohf")) uses ROHF instead, or ROKS with [rpa] xc (examples/oh-ugw-rohf.toml), semi-canonicalized per spin as for U-PDEP-RPA.

Accuracy. Smoke; treat results as about ±0.3 eV.

Limits. The quasiparticle equation is solved by a Newton root search on the self-energy, which is fragile near \( \Sigma_c \) poles. Runs report whether each root and each eigenvalue-self-consistency loop converged, and warn when one did not; check those flags.

BSE-TDA

What it is. Bethe–Salpeter excitation energies in the Tamm–Dancoff approximation on top of G0W0@HF quasiparticle energies, closed shell.

Run it. method.kind = "bse-tda" (examples/water-bse-tda.toml, and a set of *-bse-tda-augdz.toml examples for small organics); Python ferric.run_bse_tda.

Accuracy. Smoke. Given the same quasiparticle energies, the lowest five singlet excitation energies match an independent numpy BSE-TDA (PySCF density-fitted integrals, static RPA W) to 1.8e-10 Ha and their oscillator strengths to 2.6e-9, for H2O at cc-pVDZ and aug-cc-pVDZ and NH3 and CH2O at cc-pVDZ. The quasiparticle energies come from the internal G0W0@HF, which matches PySCF only at [rpa] n_quad = 100 and trunc_thresh = 0; the defaults are coarser. The core and high-virtual quasiparticle energies are ill-conditioned, which moves the lowest five excitations by at most 2.2e-7 Ha.

TDDFT and TDA

What it is. Linear-response excitations in the Tamm–Dancoff approximation (TDA, which is CIS for an HF reference) and the full Casida equations, closed shell.

Run it. method.kind = "tda" or "tddft" with [tddft] n_roots and optionally xc (examples/water-tda.toml, examples/water-tddft-pbe.toml); Python ferric.run_tddft(mol, bs, aux, functional=..., method="tda") or method="casida".

Scope. Closed-shell references, singlet excitations. With a DFT reference the \( (ia|f_{xc}|jb) \) XC-kernel term is included; with no functional the result is CIS/TDHF. Meta-GGA, VV10 and range-separated functionals are refused. Grade: Proven (narrow, closed shell): water, formaldehyde and NH3 at 6-31G and aug-cc-pVDZ with HF, LDA, PBE and B3LYP match PySCF TDA/TDDFT to at most 6.5e-4 eV, test bar 1e-3 eV (ferric-tddft/tests/validation_tddft.rs).

A separate, library-only TDA-DFT in ferric-gw/src/tddft.rs uses the same kernel (ferric_dft::lr_kernel) and is pinned against PySCF (ferric-gw/tests/tda_dft_vs_pyscf.rs). The user-facing TDA reproduces it to 1e-8 Ha for HF, PBE and B3LYP.

Polarizabilities and dispersion coefficients

What exists. Static molecular and atom-partitioned polarizabilities, Casimir–Polder \( C_6 \) coefficients from per-atom dynamic polarizabilities \( \alpha^A(i\omega) \), and many-body dispersion (MBD). Three sources feed the \( C_6 \) contraction ([rpa] c6_source):

• ts: the Tkatchenko–Scheffler single-pole model (default).
• mbd: self-consistent dipole screening of the TS polarizabilities. This is the full-range screening of Tkatchenko et al. 2012 (Gaussian-damped dipole tensor, no Fermi range separation), not the range-separated screening of MBD@rsSCS.
• pdep: dynamic PDEP-RPA polarizabilities (examples/water-c6-pdep.toml, examples/argon-c6-rpa-pbe.toml).

MBD@rsSCS dispersion energy. ferric.mbd_rsscs_energy(mol, volume_ratios, functional=... | beta=...) (Rust: ferric_rpa::dispersion::mbd_rsscs_energy) computes the MBD@rsSCS energy of Ambrosetti et al. 2014 for a finite molecule, standalone (no SCF). Inputs are per-atom Hirshfeld volume ratios, which scale the free-atom α, C6 and R_vdW (Z = 1–54). The polarizabilities are screened with the short-range part (1 − f) of the Gaussian-damped dipole tensor on libMBD's 15-point imaginary-frequency grid; the energy couples the screened oscillators through the long-range part f of the bare dipole tensor, with the Fermi function f(R) = 1/(1 + exp(−6(R/(β(R_A + R_B)) − 1))). β is functional dependent: PBE 0.83, PBE0 and HSE06 0.85; any other functional needs an explicit beta. It returns the energy and the screened α₀, C6, R_vdW and ω per atom, and raises on a polarization catastrophe (a non-positive coupled-oscillator eigenvalue). No periodic systems. As a dispersion correction on a Kohn–Sham SCF it is [dft] dispersion = "mbd" / run_dft(dispersion="mbd"), which take the volume ratios from the converged density (reported as DftResult.volume_ratios) and add the energy and its analytic gradient (see SCF and DFT).

The argon-c6-rpa-pbe.toml header records C6(Ar–Ar) = 56.4 a.u. at RPA@PBE/aug-cc-pVTZ against the DOSD value 64.3 (−12%). Which of TS and PDEP-RPA gives better molecular \( C_6 \) is not established.

Use an augmented basis for any polarizability or \( C_6 \): without diffuse functions the dipole response is badly underestimated.

TDHF/RPAx \( C_6 \) is a measured negative. \( C_6 \) built on the RPAx@PBE kernel stays about 63% low regardless of the gap. Its static polarizability is not established either: at a physical scissor (0.36 Ha) water/cc-pVDZ gives an isotropic α of 5.20 a.u. against the DOSD 9.64 a.u. (−46%). The 9.24 a.u. quoted in the example's header comes from scissor = 0, where the tensor has a negative diagonal component, and that setting is refused. method.kind = "tdhf-static-polarizability" computes static α only. It needs a KS reference ([rpa] xc), and at the default [gw] scissor = 0 it can hit an excitonic instability, which is reported as an error rather than a negative α; examples/water-tdhf-static-alpha.toml hits it as shipped, so set scissor to about 0.3–0.4 Ha.

Cite

PDEP: Wilson, Gygi & Galli 2008. RI-RPA quadrature: Eshuis, Yarkony & Furche 2010; minimax grids: Kaltak, Klimeš & Kresse 2014. GW: Hedin 1965; GW100: van Setten et al. 2015. TDDFT review: Dreuw & Head-Gordon 2005. TS: Tkatchenko & Scheffler 2009; MBD: Tkatchenko et al. 2012; MBD@rsSCS: Ambrosetti et al. 2014. Full entries in References.