OpenRDM An open-source library for reduced-density matrix-based analysis and computation
Theory

# Background

Here, we use the conventional notation of multireference (MR) methods when labeling the molecular orbitals (MOs), $$\lbrace \psi\rbrace$$: the indices $$i$$, $$j$$, $$k$$, and $$l$$ denote inactive (doubly occupied) orbitals; $$t$$, $$u$$, $$v$$, and $$w$$ represent active orbitals; and $$p$$, $$q$$, $$r$$, and $$s$$ indicate general orbitals. Einstein's summation rule is implied in all expressions.

We begin by defining the non-relativistic Born-Oppenheimer (BO) electronic Hamiltonian

$$\tag{1}\label{EQ:HAMILTONIAN} \hat{H} = h^p_q \hat{E}^p_q + \frac{1}{2} \nu^{pq}_{rs} \hat{e}^{pq}_{rs},$$

in terms of one- and two-particle excitation operators which are defined as

\begin{gather} \hat{E}^{p}_{q} = \hat{a}^\dagger_{p_\sigma} \hat{a}_{q_\sigma} \tag{2a}\label{EQ:Epq} \\ \hat{e}^{p r}_{q s} = \hat{E}^p_q \hat{E}^r_s - \delta^q_r \hat{E}^p_s = \hat{a}^\dagger_{p_\sigma} \hat{a}^\dagger_{r_\tau} \hat{a}_{s_\tau} \hat{a}_{q_\sigma} \tag{2b}\label{EQ:epqrs} \end{gather}

In Eqs. $$\eqref{EQ:Epq}$$ and $$\eqref{EQ:epqrs}$$, the $$\hat{a}^\dagger$$ and $$\hat{a}$$ represent second-quantized creation and annihilation operators, respectively, and the Greek labels run over $$\alpha$$ and $$\beta$$ spins (the sum over which is implied). The symbol $$h^p_q = \left<\psi_p|\hat{h}|\psi_q\right>$$ represents the sum of the electron kinetic energy and electron-nucleus potential energy integrals, and $$\nu^{pq}_{rs} = \left<\psi_p \psi_q|\psi_r\psi_s\right>$$ is an element of the two-electron repulsion integral tensor. Because the electronic Hamiltonian includes up to only pair-wise interactions, the ground-state energy of a many-electron system can be expressed as an exact linear functional of the the one-electron reduced-density matrix (1-RDM) and two-electron reduced-density matrix (2-RDM)

$$\tag{3}\label{EQ:Eel} E = {}^1D^p_q h^p_q + \frac{1}{2} {}^2D^{pq}_{rs} \nu^{pq}_{rs}.$$

Here, the 1-RDM and the 2-RDM are represented in their spin-free forms, with elements defined as

$$\tag{4a} {}^1D^p_q = {}^1D^{p_\sigma}_{q_\sigma} = \left<\Psi|\hat{a}^\dagger_{p_\sigma} \hat{a}_{q_\sigma}|\Psi\right> \label{EQ:1RDM}$$

and

$$\tag{4b} {}^2D^{pq}_{rs} = {}^2D^{p_\sigma q_\tau}_{r_\sigma s_\tau} = \left<\Psi|\hat{a}^\dagger_{p_\sigma} \hat{a}^\dagger_{q_\tau} \hat{a}_{s_\tau} \hat{a}_{r_\sigma}|\Psi\right> \label{EQ:2RDM},$$

respectively. Again, the summation over the spin labels in Eqs. $$\eqref{EQ:1RDM}$$ and $$\eqref{EQ:2RDM}$$ is implied.

# Multiconfiguration Pair-Density Functional Theory

The multiconfiguration pair-density functional theory (MC-DPFT) energy expression can be written as

$$\tag{5}\label{EQ:EMC-PDFT} E_{\text{MC-PDFT}} = 2h^i_i + h^t_u {}^1D^t_u + E_\text{H} + E_\text{xc}\left[\rho,\Pi,|\nabla\rho|,|\nabla\Pi|\right],$$

where the Hartree energy, $$E_\text{H}$$, is

$$\tag{6}\label{EQ:EHARTREE} E_\text{H} = 2 \nu^{ij}_{ij} + 2\nu^{ti}_{ui} {}^1D^t_u + \frac{1}{2} \nu^{tv}_{uw} {}^1D^{t}_{u} {}^1D^{v}_{w}$$

The total electronic density and its gradient are defined by the 1-RDM as

$$\label{EQ:RHO}\tag{7} \rho(\mathbf{r}) = {}^1D^p_q\ \psi^*_p(\mathbf{r}) \psi_q(\mathbf{r}),$$

and

$$\label{EQ:DRHO}\tag{8} \nabla\rho(\mathbf{r}) = {}^1D^p_q \left[ \nabla\psi^*_p(\mathbf{r}) \psi_q(\mathbf{r}) + \psi^*_p(\mathbf{r}) \nabla\psi_q(\mathbf{r}) \right],$$

respectively. The on-top pair density (OTPD) and its gradient can similarly be defined in terms of the 2-RDM as

$$\label{EQ:PI}\tag{9} \Pi(\mathbf{r}) = {}^2D^{pq}_{rs}\ \psi^*_p(\mathbf{r}) \psi^*_q(\mathbf{r}) \psi_r(\mathbf{r}) \psi_s(\mathbf{r}),$$

and

\begin{eqnarray} \label{EQ:DPI}\tag{10} \nabla\Pi(\mathbf{r}) = {}^2D^{pq}_{rs} &[& \nabla\psi^*_p(\mathbf{r}) \psi^*_q(\mathbf{r}) \psi_r(\mathbf{r}) \psi_s(\mathbf{r}) \nonumber \\\ &+& \psi^*_p(\mathbf{r}) \nabla\psi^*_q(\mathbf{r}) \psi_r(\mathbf{r}) \psi_s(\mathbf{r}) \nonumber \\\ &+& \psi^*_p(\mathbf{r}) \psi^*_q(\mathbf{r}) \nabla\psi_r(\mathbf{r}) \psi_s(\mathbf{r}) \nonumber \\\ &+& \psi^*_p(\mathbf{r}) \psi^*_q(\mathbf{r}) \psi_r(\mathbf{r}) \nabla\psi_s(\mathbf{r}) ~], \end{eqnarray}

respectively. Here, the 1- and 2-RDMs are obtained from an MR computation.

At this stage, we must identify a suitable OTPD functional for use in MC-PDFT. The simplest class of functionals can be derived from existing approximate exchange-correlation (XC) functionals employed within Kohn-Sham DFT by first recognizing that, for a density derived from a single Slater determinant, the spin magnetization can be expressed exactly in terms of the OTPD and the total density. [1, 2] More specifically, the spin polarization factor, $$\zeta(\mathbf{r}) = m(\mathbf{r})/\rho(\mathbf{r})$$, can be expressed as

$$\label{EQ:ZETR}\tag{11} \zeta(\mathbf{r}) = \sqrt{1-R(\mathbf{r})},$$

where

$$\label{EQ:OTR}\tag{12} R(\mathbf{r}) = \frac{4~ \Pi(\mathbf{r})}{\rho^2(\mathbf{r})}.$$

The basic assumption underlying the "translated" (t) OTPD functionals proposed in Ref. [3] is that the spin polarization factor can be similarly defined for a density and OTPD obtained from a MR method, as

\begin{eqnarray} \label{EQ:TR}\tag{13} \zeta_{\text{tr}}(\mathbf{r}) =& \begin{cases} \sqrt{1-R(\mathbf{r})} ,& \quad R(\mathbf{r}) \leq 1 \\\ 0 ,& \quad R(\mathbf{r}) > 1 \end{cases} \end{eqnarray}

where the second case accounts for the fact that the argument of the square root can become negative for $$\rho(\mathbf{r})$$ and $$\Pi(\mathbf{r})$$ that are not derived from a single-configuration wave function. The translated OTPD functional is then defined as

$$\label{EQ:TRE}\tag{14} E_{\text{OTPD}}\left[\rho(\mathbf{r}), \Pi(\mathbf{r}),| \nabla\rho(\mathbf{r})| \right] \equiv E_{\text{xc}}[\tilde{\rho}_\alpha(\mathbf{r}), \tilde{\rho}_\beta(\mathbf{r}), |\nabla\tilde{\rho}_\alpha(\mathbf{r})|, |\nabla\tilde{\rho}_\beta(\mathbf{r})|],$$

where the tilde refers to translated densities and their gradients, given by [3,4]

$$\label{EQ:TRHO}\tag{15} \tilde{\rho}_\sigma(\mathbf{r}) = \frac{\rho(\mathbf{r})}{2} \left(1 + c_\sigma \zeta_{\text{tr}}(\mathbf{r})\right),$$

and

$$\tag{16} \nabla\tilde{\rho}_\sigma(\mathbf{r}) = \frac{\nabla\rho(\mathbf{r})}{2} \left(1 + c_\sigma \zeta_{\text{tr}}(\mathbf{r})\right),$$

respectively. Here, $$c_\sigma$$ = 1~(-1) when $$\sigma = \alpha$$ ( $$\beta$$).

It is important to note that, in deriving the translated OTPD functional expression in Eq. $$\eqref{EQ:TRE}$$, no dependence on $$\nabla\Pi(\mathbf{r})$$ is assumed. A scheme in which the OTPD functional depends explicitly upon $$\nabla\Pi(\mathbf{r})$$ has also been proposed.[5] The corresponding "fully-translated" (ft) functionals are defined as

$$\label{EQ:FTRE}\tag{17} E_{\text{OTPD}}\left[\rho(\mathbf{r}), \Pi(\mathbf{r}), |\nabla\rho(\mathbf{r})|,|\nabla\Pi(\mathbf{r}) \right|] \equiv E_{\text{xc}}[\tilde{\rho}_\alpha(\mathbf{r}), \tilde{\rho}_\beta(\mathbf{r}), |\nabla\tilde{\rho}_\alpha(\mathbf{r})|, |\nabla\tilde{\rho}_\beta(\mathbf{r})|]$$