13.2 Gradients Theory

All details on the theory and results are published in [300, 301]. The full-core and frozen-core treatments can be handled using the same general formalism described below. The RI-RPA energy is a function of the MO coefficients $\mathbf{C}$ and the Lagrange multipliers $\boldsymbol{\epsilon}$ and depends parametrically (i) on the interacting Hamiltonian $\hat{H}$, (ii) on the AO basis functions and the auxiliary basis functions. All parameters may be gathered in a supervector $\mathbf{X}$ and thus

\[\begin{equation} E^{\text{RIRPA}} \equiv E^{\text{RIRPA}}(\mathbf{C},\boldsymbol{\epsilon}|\mathbf{X}) \;\text{.} \end{equation}\](13.9)

$\mathbf{C}$ and $\boldsymbol{\epsilon}$ in turn depend parametrically on $\mathbf{X}$, the exchange-correlation matrix $\mathbf{V}^{\text{XC}}$, and the overlap matrix $\mathbf{S}$ through the KS equations and the orbital orthonormality constraint. First-order properties may be defined in a rigorous and general fashion as total derivatives of the energy with respect to a “perturbation” parameter $\xi$. However, the RI-RPA energy is not directly differentiated in our method. Instead, we define the RI-RPA energy Lagrangian

\[\begin{align} &\quad L^{\text{RIRPA}}(\mathcal{C},\mathcal{E},\mathcal{D}^\Delta,\mathcal{W}|\mathbf{X},\mathbf{V}^{\text{XC}},\mathbf{S}) \\ &= E^{\text{RIRPA}}(\mathcal{C},\mathcal{E}|\mathbf{X}) + \sum_\sigma \left(\left\langle \mathcal{D}_\sigma^\Delta(\mathcal{C}_\sigma^{\text{T}} \mathbf{F}_\sigma \mathcal{C}_\sigma - \mathcal{E}_\sigma)\right\rangle - \left\langle \mathcal{W}_\sigma (\mathcal{C}_\sigma^{\text{T}} \mathbf{S} \mathcal{C}_\sigma - \mathbf{1}) \right\rangle \right) \;\text{.} \end{align}\](13.10)

$\mathcal{C}$, $\mathcal{E}$, $\mathcal{D}^\Delta$, and $\mathcal{W}$ are independent variables. $L^{\text{RIRPA}}$ is required to be stationary with respect to $\mathcal{C}$, $\mathcal{E}$, $\mathcal{D}^\Delta$, and $\mathcal{W}$. $\mathcal{D}^\Delta$ and $\mathcal{W}$ act as Lagrange multipliers enforcing that $\mathcal{C}$ and $\mathcal{E}$ satisfy the KS equations and the orbital orthonormality constraint,

\[\begin{align} \left(\frac{\partial L^{\text{RIRPA}}}{\partial \mathcal{D}_\sigma^\Delta} \right)_{\text{stat}} &= \mathcal{C}_\sigma^{\text{T}} \mathbf{F}_\sigma \mathcal{C}_\sigma - \mathcal{E}_\sigma = \mathbf{0} \;\text{,} \\ \left(\frac{\partial L^{\text{RIRPA}}}{\partial \mathcal{W}_\sigma} \right)_{\text{stat}} &= \mathcal{C}_\sigma^{\text{T}} \mathbf{S} \mathcal{C}_\sigma - \mathbf{1} = \mathbf{0} \;\text{.} \end{align}\](13.11–13.12)

$\mathcal{D}^\Delta$ and $\mathcal{W}$ are determined by the remaining stationarity conditions,

\[\begin{equation} \left( \frac{\partial L^{\text{RIRPA}}}{\partial \mathcal{E}} \right)_{\text{stat}} = \mathbf{0} \;\text{,} \end{equation}\](13.13)

and

\[\begin{equation} \left( \frac{\partial L^{\text{RIRPA}}}{\partial \mathcal{C}} \right)_{\text{stat}} = \mathbf{0} \;\text{.} \end{equation}\](13.14)

It turns out from eqs 13.13 and 13.14 that the determination of $\mathcal{D}^\Delta$ and $\mathcal{W}$ requires the solution of a single Coupled-Perturbed KS equation. Complete expressions for $\mathcal{D}^\Delta$ and $\mathcal{W}$ are given in [300]. At the stationary point “$\text{stat}=(\mathcal{C}=\mathbf{C},\mathcal{E}=\boldsymbol{\epsilon},\mathcal{D}^\Delta=\mathbf{D}^\Delta,\mathcal{W}=\mathbf{W})$”, first-order RI-RPA properties are thus efficiently obtained from

\[\begin{align} \frac{dE^{\text{RIRPA}}(\mathbf{C},\boldsymbol{\epsilon}|\mathbf{X})}{d\xi} &= \left\langle \left(\frac{\partial L^{\text{RIRPA}}}{\partial \mathbf{X}}\right)_{\text{stat}} \frac{d\mathbf{X}}{d\xi}\right\rangle + \left\langle \left(\frac{\partial L^{\text{RIRPA}}}{\partial \mathbf{V}^{\text{XC}}}\right)_{\text{stat}} \left(\frac{\partial\mathbf{V}^{\text{XC}}}{\partial\xi}\right)_{\text{stat}}\right\rangle \\ &\quad + \left\langle\left(\frac{\partial L^{\text{RIRPA}}}{\partial \mathbf{S}}\right)_{\text{stat}} \frac{d\mathbf{S}}{d\xi} \right\rangle \;\text{.} \end{align}\](13.15)

Finally, the RPA energy gradients may be explicitly expanded as follows:

\[\begin{align} \frac{d E^{\text{RIRPA}}(\mathbf{C},\boldsymbol{\epsilon}|\mathbf{X})}{d\xi} &= \left\langle \mathbf{D}^{\text{RIRPA}} \frac{d\mathbf{h}}{d\xi} \right\rangle + \left\langle \boldsymbol{\Gamma}^{(4)} \frac{d\boldsymbol{\Pi}^{(4)}}{d\xi} \right\rangle + \left\langle \mathbf{D}^{\Delta} \left(\frac{\partial \mathbf{V}^\text{XC}[\mathbf{D}]}{\partial\xi}\right)_\text{stat} \right\rangle \\ &\quad + \left\langle \boldsymbol{\Gamma}^{(3)} \frac{d\boldsymbol{\Pi}^{(3)}}{d\xi} \right\rangle + \left\langle \boldsymbol{\Gamma}^{(2)} \frac{d\boldsymbol{\Pi}^{(2)}}{d\xi} \right\rangle - \left\langle \mathbf{W} \frac{d\mathbf{S}}{d\xi} \right\rangle \;\text{.} \end{align}\](13.16)

where $\mathbf{D}^{\text{RIRPA}}$ is the KS ground state one-particle density matrix $\mathbf{D}$ plus the RI-RPA difference density matrix $\mathbf{D}^{\Delta}$ which corrects for correlation and orbital relaxation effects. $\mathbf{h}$ is the one-electron Hamiltonian; $\boldsymbol{\Pi}^{(2/3/4)}$ are 2-, 3-, and 4-centre electron repulsion integrals and the $\boldsymbol{\Gamma}^{(2/3/4)}$ are the corresponding 2-, 3-, and 4-index relaxed 2-particle density matrices; $\mathbf{W}$ may be interpreted as the energy-weighted total spin one-particle density matrix.

This result illustrates the key advantage of the Lagrangian method: Total RI-RPA energy derivatives featuring a complicated implicit dependence on the parameter $\mathbf{X}$ through the variables $\mathbf{C}$ and $\boldsymbol{\epsilon}$ are replaced by partial derivatives of the RI-RPA Lagrangian, whose computation is straightforward once the stationary point of the Lagrangian has been fully determined.

Gradients Prerequisites

Geometry optimizations and first order molecular property calculations can be executed by adding the keyword rpagrad to the $rirpa section in the control file. RPA gradients also require

an auxiliary basis defined in the data group $jbas for the computation of the Coulomb integrals for the Hartree-Fock energy
an auxiliary basis defined in the data group $cbas for the ERI’s in the correlation treatment.

The following gradient-specific options may be further added to the $rirpa section in the control file

drimp2 - computes gradients in the DRIMP2 limit.
niapblocks $\langle \text{integer} \rangle$ - Manual setting of the integral block size in subroutine rirhs.f; for developers.

In order to run a geometry optimization, jobex must be invoked with the level set to rirpa, and the -ri option (E.g. jobex -ri -level rirpa).

In order to run a numerical frequency calculation, NumForce must be invoked with the level set to rirpa, e.g., NumForce -d 0.02 -central -ri -level rirpa.