FCI

Given a set of \(M\) spin orbitals, the full \(N\)-electron CI wave function is written as a linear combination of all Slater determinants of all Slater determinants belonging to the Fock space \(\mathcal{F}(M,N)\): \[\begin{equation} \require{physics} \ket{\text{FCI}(\vb{c})} % = \sum_{\vb{k}}^{\mathcal{F}(M,N)} c_{\vb{k}} \ket{\vb{k}} % \thinspace . \end{equation}\] We can equally write this as \[\begin{equation} \label{eq:FCI_1_plus_C} \ket{\text{FCI}(\vb{c})} % = \qty(1 + \hat{C}) \ket{\text{HF}} % \thinspace , \end{equation}\] in which \(\hat{C}\) is a linear excitation operator: \[\begin{equation} \hat{C} % = \sum_I^{N} \sum_{A=N+1}^M % c_I^A \hat{E}_{AI} % + \sum_{IJ}^{N} \sum_{AB=N+1}^M % c_{IJ}^{AB} \hat{E}_{JB} \hat{E}_{AI} % + \cdots % \thinspace . \end{equation}\]

The FCI wave function is the electronic state that is an eigenfunction of the Hamiltonian: \[\begin{equation} \hat{\mathcal{H}} \ket{\text{FCI}} = E \ket{\text{FCI}} % \thinspace . \end{equation}\] We can, in theory, solve this equation by expanding in the ONV basis of the Fock space, leading to the matrix eigenvalue problem \[\begin{equation} \vb{H} \vb{C} = E \vb{C} \thinspace , \end{equation}\] in which \(\vb{H}\) is the Hamiltonian matrix in the ONV basis: \[\begin{equation} \vb{H}_{\vb{k} \vb{m}} % = \matrixel{ \vb{k} }{ \hat{\mathcal{H}} }{ \vb{m} } % \thinspace , \end{equation}\] and \(\vb{C}\) is the coefficient vector in the ONV basis: \[\begin{equation} \vb{C}_{\vb{m}} = \braket{\vb{m}}{\text{FCI}} % \thinspace , \end{equation}\] in which \(\vb{k}\) and \(\vb{m}\) are labels for the \((\dim{\mathcal{F}(M,N)})\)-dimensional ONV basis vector.

In this section, we will derive/list FCI Hamiltonian matrix elements, which are all of the form \[\begin{equation} \mathcal{H}_{I_\alpha I_\beta, J_\alpha J_\beta} % = \matrixel{ I_\alpha I_\beta % }{ \hat{\mathcal{H}} }{ % J_\alpha J_\beta % } % \tag{\ref{eq:FCI_Hamiltonian_matrix_elements}} % \thinspace , \end{equation}\] making use of the Slater-Condon rules. Filling in the expression for the molecular Hamiltonian, making use of the effective one-electron integrals from equation \(\eqref{eq:effective_one_electron_integrals}\) and finally splitting the singlet one-electron excitation operators in their \(\alpha\) and \(\beta\) parts (cfr. equation \(\eqref{eq:E_pq_alpha_beta}\)), we get \[\begin{align} \mathcal{H}_{I_\alpha I_\beta, J_\alpha J_\beta} % &= \sum_{pq}^K % k_{pq} % \matrixel{ % I_\alpha I_\beta % }{ \hat{E}_{pq} }{ % J_\alpha J_\beta % } % + \frac{1}{2} \sum_{pqrs}^K g_{pqrs} % \matrixel{ % I_\alpha I_\beta % }{ \hat{E}_{pq} \hat{E}_{rs} }{ % J_\alpha J_\beta % } \\ &= \mathcal{H}^{(1)}_{I_\alpha I_\beta, J_\alpha J_\beta} % + \mathcal{H}^{(2)}_{I_\alpha I_\beta, J_\alpha J_\beta} \\ &= \qty( % \mathcal{H}^{\alpha}_{I_\alpha I_\beta, J_\alpha J_\beta} % + \mathcal{H}^{\beta}_{I_\alpha I_\beta, J_\alpha J_\beta} % ) % \qty( % \mathcal{H}^{\alpha \alpha}_{I_\alpha I_\beta, J_\alpha J_\beta} % + \mathcal{H}^{\alpha \beta}_{I_\alpha I_\beta, J_\alpha J_\beta} % + \mathcal{H}^{\beta \beta}_{I_\alpha I_\beta, J_\alpha J_\beta} % ) % \thinspace . \end{align}\]

For the diagonal (i.e. \(I = J\)) one-electron contributions, we have: \[\begin{align} &\mathcal{H}^{\alpha}_{I_\alpha I_\beta, I_\alpha I_\beta} % = \sum_{p \in I_\alpha}^K k_{pp} \\ &\mathcal{H}^{\beta}_{I_\alpha I_\beta, I_\alpha I_\beta} % = \sum_{p \in I_\beta}^K k_{pp} \thinspace , \end{align}\] and for the diagonal two-electron contributions, we have: \[\begin{align} & \mathcal{H}^{\alpha \alpha}_{I_\alpha I_\beta, I_\alpha I_\beta} % = \frac{1}{2} \sum_{p,q \in I_\alpha}^K g_{ppqq} % + \frac{1}{2} \sum_{p \in I_\alpha}^K % \sum_{q \not\in I_\alpha}^K g_{pqqp} \\ % & \mathcal{H}^{\alpha \beta}_{I_\alpha I_\beta, I_\alpha I_\beta} % = \sum_{p \in I_\alpha}^K \sum_{q \in I_\beta}^K g_{ppqq} \\ % & \mathcal{H}^{\beta \beta}_{I_\alpha I_\beta, I_\alpha I_\beta} % = \frac{1}{2} \sum_{p,q \in I_\beta}^K g_{ppqq} % + \frac{1}{2} \sum_{p \in I_\beta}^K % \sum_{q \not\in I_\beta}^K g_{pqqp} % \thinspace . \end{align}\]

If we have \(\ket{I_\alpha}\) and \(\ket{J_\alpha}\) that differ in one electron excitation, for example in the orbitals \(w \alpha\) and \(x \alpha\), i.e. \[\begin{align} &\ket{J_\alpha} % = \ket{ k_1, \ldots, 0_{w \alpha}, \ldots, 1_{x \alpha}, \ldots, k_{K \alpha} } \\ &\ket{I_\alpha} = \ket{ k_1, \ldots, 1_{w \alpha}, \ldots, 0_{x \alpha}, \ldots, k_{K \alpha} } \end{align}\] and similarly for \(\beta\) string, we have for the one-electron part: \[\begin{align} & \mathcal{H}^{\alpha}_{I_\alpha I_\beta, J_\alpha J_\beta} % = \delta_{I_\beta J_\beta} % \Gamma^{I_\alpha}_{w \alpha} % \Gamma^{J_\alpha}_{x \alpha} % \thinspace k_{wx} \\ & \mathcal{H}^{\beta}_{I_\alpha I_\beta, J_\alpha J_\beta} % = \delta_{I_\alpha J_\alpha} % \Gamma^{I_\beta}_{w \beta} % \Gamma^{J_\beta}_{x \beta} % \thinspace k_{wx} \thinspace , \end{align}\] and for the two-electron part: \[\begin{align} & \mathcal{H}^{\alpha \alpha}_{I_\alpha I_\beta, J_\alpha, J_\beta} % = \delta_{I_\beta J_\beta} % \Gamma^{I_\alpha}_{w \alpha} % \Gamma^{J_\alpha}_{x \alpha} \qty( % \frac{1}{2} \sum_q^K g_{wqqx} % + \sum_{p \in I_\alpha}^K (g_{wxpp} - g_{wppx}) % ) \\ % & \mathcal{H}^{\alpha \beta}_{I_\alpha I_\beta, J_\alpha J_\beta} % = \Gamma^{I_\alpha}_{w \alpha} % \Gamma^{J_\alpha}_{x \alpha} % \Gamma^{I_\beta}_{w \beta} % \Gamma^{J_\beta}_{x \beta} % \thinspace g_{wxwx} \\ & \mathcal{H}^{\beta \beta}_{I_\alpha I_\beta, J_\alpha, J_\beta} % = \delta_{I_\alpha J_\alpha} % \Gamma^{I_\beta}_{w \beta} % \Gamma^{J_\beta}_{x \beta} \qty( % \frac{1}{2} \sum_q^K g_{wqqx} % + \sum_{p \in I_\beta}^K (g_{wxpp} - g_{wppx}) % ) \thinspace . \end{align}\] Recognizing the effective one-electron integrals, we can collect both the \(\alpha \alpha\) and \(\beta \beta\) one-electron part and the two-electron part for ONVs that differ in one electron excitation, resulting in \[\begin{align} & \mathcal{H}^{\alpha}_{I_\alpha I_\beta, J_\alpha J_\beta} % + \mathcal{H}^{\alpha \alpha}_{I_\alpha I_\beta, J_\alpha, J_\beta} % = \delta_{I_\beta J_\beta} % \Gamma^{I_\alpha}_{w \alpha} % \Gamma^{J_\alpha}_{x \alpha} \qty( % h_{wx} + \sum_{p \in I_\alpha}^K (g_{wxpp} - g_{wppx}) % ) \\ & \mathcal{H}^{\beta}_{I_\alpha I_\beta, J_\alpha J_\beta} % + \mathcal{H}^{\beta \beta}_{I_\alpha I_\beta, J_\alpha, J_\beta} % = \delta_{I_\alpha J_\alpha} % \Gamma^{I_\beta}_{w \beta} % \Gamma^{J_\beta}_{x \beta} \qty( % h_{wx} + \sum_{p \in I_\beta}^K (g_{wxpp} - g_{wppx}) % ) \thinspace . \end{align}\]

Let’s say we have \(\ket{I_\alpha}\) and \(\ket{J_\alpha}\) that differ in two electron excitations, for example in the orbitals \(w \alpha\), \(x \alpha\), \(y \alpha\) and \(z \alpha\), i.e. \[\begin{align} &\ket{J_\alpha} % = \ket{ k_1, \ldots, 0_{w \alpha}, \ldots, 0_{x \alpha}, \ldots, 1_{y \alpha}, \ldots, 1_{z \alpha}, \ldots, k_{K \alpha} } \\ & \ket{I_\alpha} % = \ket{ k_1, \ldots, 1_{w \alpha}, \ldots, 1_{x \alpha}, \ldots, 0_{y \alpha}, \ldots, 0_{z \alpha}, \ldots, k_{K \alpha} } % \thinspace , \end{align}\] with \(w < x\) and \(y < z\) (if not, a simple relabeling suffices), and similarly for the \(\beta\) strings. By the Slater-Condon rules, the one-electron contributions vanish, and the two-electron contributions are \[\begin{align} & \mathcal{H}^{\alpha \alpha}_{I_\alpha I_\beta, J_\alpha J_\beta} % = \delta_{I_\beta J_\beta} % \Gamma^{I_\alpha}_{w \alpha} % \Gamma^{I_\alpha}_{x \alpha} % \Gamma^{J_\alpha}_{y \alpha} % \Gamma^{J_\alpha}_{z \alpha} (g_{wyxz} - g_{wzxy}) \\ % &\mathcal{H}^{\beta \beta}_{I_\alpha I_\beta, J_\alpha J_\beta} % = \delta_{I_\alpha J_\alpha} % \Gamma^{I_\beta}_{w \beta} % \Gamma^{I_\beta}_{x \beta} % \Gamma^{J_\beta}_{y \beta} % \Gamma^{J_\beta}_{z \beta} (g_{wyxz} - g_{wzxy}) % \thinspace , \end{align}\] in which the \(\alpha\beta\)-contribution vanishes, as the disconnected \(\alpha\)- and \(\beta\)-operators are effectively one-electron operators such that their matrix elements between ONVs differing in two electron excitations vanish.

However, let’s suppose that we have some kind of reference state \(\ket{0}\): \[\begin{equation} \ket{0} = \sum_{\vb{k}} C^{(0)}_{\vb{k}} \ket{\vb{k}} \end{equation}\] that is normalized: \[\begin{align} &\braket{0} = 1 \\ &\vb{C}^{(0)} \cdot \vb{C}^{(0)} = 1 \\ &\sum_{\vb{k}} (C^{(0)}_{\vb{k}})^2 = 1 \thinspace . \end{align}\] This reference state can be the Hartree-Fock wave function (i.e. there is only one coefficient different from zero and that is \(C_{\text{HF}} = 1\)), or it could be a previous iteration of the iterative CI solver. We are thus essentially writing the FCI wave function as the normalized vector \[\begin{align} \label{eq:CI_parametrization1} & \ket{\text{FCI}} % = \frac{ % \ket{0} + \ket{\vb{c}} % }{ % \sqrt{ (\vb{C}^{(0)} + \vb{c}) \cdot (\vb{C}^{(0)} + \vb{c}) } % } \\ & \vb{C} % = \frac{ % \vb{C}^{(0)} + \vb{c} % }{ % \sqrt{( \vb{C}^{(0)} + \vb{c}) \cdot (\vb{C}^{(0)} + \vb{c})} % } \thinspace , \end{align}\] in which \(\ket{\vb{c}}\) is a small variation on \(\ket{0}\). Unfortunately, all FCI states of the form \[\begin{equation} \ket{\vb{C}_\alpha} = \ket{\vb{c}} + \alpha ( \ket{0} + \ket{\vb{c}} ) \end{equation}\] are actually independent of \(\alpha\), as \(\ket{\vb{C}_\alpha} = \ket{\text{FCI}}\) for any real \(\alpha\). We therefore call \(( \ket{0} + \vb{c} )\) the redundant component of the parametrization \(\eqref{eq:CI_parametrization1}\). Since it includes the variations \(\vb{c}\) themselves, this is a bad thing and should be avoided. Fortunately, a solution is quickly provided by introducing the projection operator \[\begin{align} &\hat{P} = 1 - \dyad{0} \\ &\vb{P} = \vb{I} - \vb{C}^{(0)} \vb{C}^{(0) \text{T}} % \thinspace , \end{align}\] such that the FCI wave function is then parametrized as \[\begin{equation} \label{eq:FCI_parametrization2} \ket{\text{FCI}} % = \frac{ % \ket{0} + \hat{P} \ket{\vb{c}} % }{ % \sqrt{ 1 + \matrixel{\vb{c}}{\hat{P}}{\vb{c}} } % } \thinspace . \end{equation}\]

We can use the previous parametrization of the FCI wave function \(\eqref{eq:FCI_parametrization2}\) to calculate the electronic energy: \[\begin{align} E(\vb{C}^{(0)} + \vb{c}) % &= \frac{ % \matrixel{0}{ \hat{\mathcal{H}} }{0} % + 2 \matrixel{\vb{c}}{ % \hat{P} \hat{\mathcal{H}} % }{0} % + \matrixel{\vb{c}}{ % \hat{P} \hat{\mathcal{H}} \hat{P} % }{\vb{c}} % }{ % 1 + \matrixel{\vb{c}}{ \hat{P} }{\vb{c}} % } \\ &= \frac{ % E(\vb{C}^{(0)}) % + 2 \vb{c}^\text{T} \vb{P} \vb{H} \vb{C}^{(0)} % + \vb{c}^\text{T} \vb{P} \vb{H} \vb{P} \vb{c} % }{ % 1 + \vb{c}^\text{T} \vb{P} \vb{c} % } \thinspace , \end{align}\] where in the last line we have used the ONV basis expansion. We can then calculate the electronic gradient as \[\begin{equation} \vb{E}^{(1)} % = 2 \vb{P} \vb{H} \vb{C}^{(0)} % = 2 (\vb{H} - E^{(0)} \vb{I}) \vb{C}^{(0)} \end{equation}\] and the the electronic Hessian as \[\begin{equation} \vb{E}^{(2)} = 2 \vb{P} (\vb{H} - E^{(0)} \vb{I}) \vb{P} % \thinspace . \end{equation}\]

Requiring the electronic gradient to vanish, we obtain the CI eigenvalue equation \[\begin{equation} \label{eq:CI_eigenvalue_equation} \vb{H} \vb{C}^{(0)} = E^{(0)} \vb{C}^{(0)} % \thinspace . \end{equation}\] This means that when we (finally, after many iterations \(\vb{C} = \vb{C}^{(0)} + \vb{c}\)) find a gradient that is zero, we rest assured that the current reference state \(\ket{0}\) given by the coefficients \(\vb{C}^{(0)}\) is an eigenfunction of the Hamiltonian.

In order to obtain CI energies, we can set up Newton’s method for the optimization of the scalar function \(E(\vb{C}^{(0)} + \vb{c})\). We will expand the energy around \(\vb{C}^{(0)}\) and terminate up to second order: \[\begin{equation} E(\vb{c}) % = E^{(0)} % + \vb{c}^\text{T} \vb{E}^{(1)} % + \frac{1}{2} \vb{c}^\text{T} \vb{E}^{(2)} \vb{c} % \thinspace , \end{equation}\] so that the Newton step is given as \[\begin{align} & \vb{E}^{(2)} \vb{c} = - \vb{E}^{(1)} \\ & \vb{P} (\vb{H} - E^{(0)} \vb{I}) \vb{P} \vb{c} % = (\vb{H} - E^{(0)} \vb{I}) \vb{C}^{(0)} % \thinspace . \end{align}\] By working out this equation, we can find the Newton step as \[\begin{equation} \vb{c} = % - \vb{C}^{(0)} % + \frac{ % (\vb{H} - E^{(0)} \vb{I})^{-1} \vb{C}^{(0)} % }{ % \vb{C}^{(0) \text{T}} (\vb{H} - E^{(0)} \vb{I})^{-1} \vb{C}^{(0)} % } \thinspace , \end{equation}\] but finding the inverse of \((\vb{H} - E^{(0)} \vb{I})\), which is generally a huge matrix, is very hard. It seems like all hope is lost, since we won’t be able to calculate \((\vb{H} - E^{(0)} \vb{I})^{-1}\), but actually we don’t have to. Davidson proposed an approximation of a quasi-Newton scheme, leading to \[\begin{equation} \vb{c} % = - (\vb{H}_0 - E^{(0)} \vb{I})^{-1} % (\vb{H} - E^{(0)} \vb{I}) % \vb{C}^{(0)} \thinspace , \end{equation}\] in which \(\vb{H}_0\) is taken to be a diagonal matrix with elements equal to those of \(\vb{H}\), such that \((\vb{H}_0 - E^{(0)} \vb{I})^{-1}\) can be trivially calculated as \(E^{(0)}\) is just the electronic energy from the initial or previous step. We can also write out the Davidson step as \[\begin{equation} \vb{c} = % - (\vb{H}_0 - E^{(0)} \vb{I})^{-1} \vb{H} \vb{C}^{(0)} % + E^{(0)} (\vb{H}_0 - E^{(0)} \vb{I})^{-1} \vb{C}^{(0)} % \thinspace , \end{equation}\] in which the computationally expensive part are matrix-vector products (matvec) of the type \(\vb{H} \vb{C}^{(0)}\).

Let us first split up the \(\boldsymbol{\sigma}\)-vector in one- and two-electron contributions: \[\begin{equation} \sigma_{I_\alpha I_\beta} % = \sigma_{I_\alpha I_\beta}^{(1)} % + \sigma_{I_\alpha I_\beta}^{(2)} \thinspace , \end{equation}\] in which \[\begin{equation} \sigma_{I_\alpha I_\beta}^{(1)} % = \sum_{J_\alpha J_\beta} \sum_{pq}^K % k_{pq} \matrixel{ % I_\alpha I_\beta % }{ \hat{E}_{pq} }{ % J_\alpha J_\beta % } C_{J_\alpha J_\beta} \end{equation}\] and \[\begin{equation} \sigma_{I_\alpha I_\beta}^{(2)} % = \frac{1}{2} \sum_{J_\alpha J_\beta} \sum_{pqrs}^K % g_{pqrs} \matrixel{ % I_\alpha I_\beta % }{ \hat{E}_{pq} \hat{E}_{rs} }{ % J_\alpha J_\beta % } C_{J_\alpha J_\beta} % \thinspace . \end{equation}\]

If we split up the one-electron excitation operator in its \(\alpha\)- and its \(\beta\)-parts, we can immediately find that \[\begin{equation} \sigma_{I_\alpha I_\beta}^{(1)} % = \sigma_{I_\alpha I_\beta}^{\alpha} % + \sigma_{I_\alpha I_\beta}^{\beta} % \thinspace , \end{equation}\] in which \[\begin{align} & \sigma_{I_\alpha I_\beta}^{\alpha} % = \sum_{J_\alpha} \sum_{pq}^K % k_{pq} \matrixel{ I_\alpha }{ \hat{E}^\alpha_{pq} }{ J_\alpha } % C_{J_\alpha I_\beta} \\ & \sigma_{I_\alpha I_\beta}^{\beta} % = \sum_{J_\beta} \sum_{pq}^K k_{pq} % \matrixel{ I_\beta }{ \hat{E}^\beta_{pq} }{ J_\beta } % C_{I_\alpha J_\beta} % \thinspace . \end{align}\]

Furthermore, we can also find \[\begin{equation} \sigma_{I_\alpha I_\beta}^{(2)} % = \sigma_{I_\alpha I_\beta}^{\alpha \alpha} % + \sigma_{I_\alpha I_\beta}^{\beta \beta} % + \sigma_{I_\alpha I_\beta}^{\alpha \beta} % \thinspace , \end{equation}\] in which \[\begin{align} & \sigma_{I_\alpha I_\beta}^{\alpha \alpha} % = \frac{1}{2} \sum_{J_\alpha} % \sum_{pqrs}^K g_{pqrs} % \matrixel{ I_\alpha }{ % \hat{E}^{\alpha}_{pq} \hat{E}^{\alpha}_{rs} % }{ J_\alpha } % C_{J_\alpha I_\beta} \\ % & \sigma_{I_\alpha I_\beta}^{\beta \beta} % = \frac{1}{2} \sum_{J_\beta} % \sum_{pqrs}^K g_{pqrs} % \matrixel{ I_\beta }{ % \hat{E}^{\beta}_{pq} \hat{E}^{\beta}_{rs} % }{ J_\beta } % C_{I_\alpha J_\beta} \\ % & \sigma_{I_\alpha I_\beta}^{\alpha \beta} % = \sum_{J_\alpha J_\beta} % \sum_{pqrs}^K g_{pqrs} % \matrixel{ I_\alpha }{ % \hat{E}^\alpha_{pq} % }{ J_\alpha } % \matrixel{ I_\beta }{ % \hat{E}^\beta_{rs} % }{ J_\beta } % C_{J_\alpha J_\beta} % \thinspace . \end{align}\]

For the diagonal contributions (i.e. \(p=q\)) to the 1-DMs, we have \[\begin{align} & D^{\alpha \alpha}_{pp} % = \sum_{I \alpha}^{\dim \mathcal{F}_\alpha} % \sum_{I \beta}^{\dim \mathcal{F}_\beta} % |c_{I_\alpha I_\beta}|^2 \\ % & D^{\beta \beta}_{pp} % = \sum_{I \alpha}^{\dim \mathcal{F}_\alpha} % \sum_{I \beta}^{\dim \mathcal{F}_\beta} % |c_{I_\alpha I_\beta}|^2 % \thinspace , \end{align}\] in which it is understood that for \(D^{\alpha \alpha}_{pp}\), \(p\) must be occupied in \(I_\alpha\) (\(p \in I_\alpha\)), and similarly for \(\beta\). For the off-diagonal elements (i.e. \(p \neq q\)), we have \[\begin{align} & D^{\alpha \alpha}_{pq} % = \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % \Gamma^{I_\alpha}_{p \alpha} % \Gamma^{J_\alpha}_{q \alpha} \qty( % \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % c^*_{I_\alpha I_\beta} % c_{J_\alpha I_\beta} % ) \\ % & D^{\beta \beta}_{pq} % = \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % \Gamma^{I_\beta}_{p \beta} % \Gamma^{J_\beta}_{q \beta} \qty( % \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % c^*_{I_\alpha I_\beta} % c_{I_\alpha J_\beta} % ) \thinspace , \end{align}\] in which, for \(D^{\alpha \alpha}_{pq}\), it is understood that \(p \in I_\alpha\) and \(q \not\in I_\alpha\) and \(\ket{J_\alpha} = \hat{a}^\dagger_{q \alpha} \hat{a}_{p \alpha} \ket{I_\alpha}\), and similarly for \(\beta\).

For the 2-DMs, we have the following formulas: \[\begin{align} & d^{\alpha \alpha \alpha \alpha}_{ppqq} % = \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % |c_{I_\alpha I_\beta}|^2 \\ & d^{\beta \beta \beta \beta}_{ppqq} % = \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % |c_{I_\alpha I_\beta}|^2 % \thinspace , \end{align}\] in which we understand that we must have \(p \in I_\alpha\) and \(q \in I_\alpha\), and similarly for \(\beta\). We also have \[\begin{align} & d^{\alpha \alpha \alpha \alpha}_{ppqr} % = \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % \Gamma^{I_\alpha}_{q \alpha} % \Gamma^{J_\alpha}_{r \alpha} \qty( % \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % c^*_{I_\alpha I_\beta} c_{J_\alpha I_\beta} % ) \\ & d^{\beta \beta \beta \beta}_{ppqr} % = \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % \Gamma^{I_\beta}_{q \beta} % \Gamma^{J_\beta}_{r \beta} \qty( % \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % c^*_{I_\alpha I_\beta} c_{I_\alpha J_\beta} % ) \thinspace , \end{align}\] in which it is understood that \(p \in I_\alpha\), \(q \in I_\alpha\) and \(r \not\in I_\alpha\), and that the string that couples to \(I_\alpha\) is \(\ket{J_\alpha} = \hat{a}^\dagger_{r \alpha} \hat{a}_{q \alpha} \ket{I_\alpha}\), and of course similarly for \(\beta\). In the case where \(p \neq q \neq r \neq s\), we have \[\begin{align} & d^{\alpha \alpha \alpha \alpha}_{pqrs} % = \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % \Gamma^{J_\alpha}_{q \alpha} % \Gamma^{J_\alpha}_{s \alpha} % \Gamma^{I_\alpha}_{p \alpha} % \Gamma^{I_\alpha}_{r \alpha} \qty( % \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % c^*_{I_\alpha I_\beta} c_{J_\alpha I_\beta} % ) \\ % & d^{\beta \beta \beta \beta}_{pqrs} = % \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % \Gamma^{J_\beta}_{q \beta} % \Gamma^{J_\beta}_{s \beta} % \Gamma^{I_\beta}_{p \beta} % \Gamma^{I_\beta}_{r \beta} \qty( % \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % c^*_{I_\alpha I_\beta} c_{I_\alpha J_\beta} % ) \thinspace , \end{align}\] in which we understand that the coupling string is \(\ket{J_\alpha} = \hat{a}^\dagger_{q \alpha} \hat{a}^\dagger_{s \alpha} \hat{a}_{r \alpha} \hat{a}_{p \alpha} \ket{I_\alpha}\), such that \(p \in I_\alpha\), \(r \in I_\alpha\), \(s \not\in I_\alpha\) and \(q \not\in I_\alpha\), and similarly for \(\beta\).

For the 2-DMs with mixed spin, we find \[\begin{equation} d^{\alpha \alpha \beta \beta}_{pppp} % = \sum_{I_\alpha}^{\dim \mathcal{F}_\alpha} % \sum_{I_\beta}^{\dim \mathcal{F}_\beta} % |c_{I_\alpha I_\beta}|^2 % \thinspace , \end{equation}\] in which \(p \in I_\alpha\) and \(p \in I_\beta\).