AP1roG PSEs
AP1roG being a special case of APIG leads to simplifications to equations \(\eqref{eq:APIG_overlap_paired}\), \(\eqref{eq:APIG_overlap_single_pair_excitation}\), and \(\eqref{eq:APIG_overlap_double_pair_excitation}\): \[\begin{equation} \require{physics} \braket{\Phi_0}{\text{AP1roG}} = 1 \thinspace , \end{equation}\] \[\begin{equation} \braket{\Phi_i^a}{\text{AP1roG}} = G^a_i \thinspace , \end{equation}\] and \[\begin{equation} \braket{\Phi_{ij}^{ab}}{\text{AP1roG}} = G^a_i G^b_j + G^a_j G^b_i \thinspace . \end{equation}\]
Using these expressions in \(\eqref{eq:PSE_APIG}\), we get: \[\begin{equation} \label{eq:pse_ap1rog} \begin{split} & g_{ai} - g_{ia} \qty(G_i^a)^2 + \qty[ 2 (h_{aa} - h_{ii}) + (g_{aa} - g_{ii}) + \sum_{j \neq i}^{N_P} (\gamma_{aj} + \gamma_{ja} - \gamma_{ij} - \gamma_{ji} ) ] G_i^a \\ & + \sum_{j \neq i}^{N_P} \Big(g_{ji} - g_{ja} G_i^a \Big) G_j^a + \sum_{b = N_P + 1, b \neq a}^K \Big( g_{ab} - g_{ib} G_i^a \Big) G_i^b + \sum_{j \neq i}^{N_P} \sum_{b = N_P + 1, b \neq a}^K g_{jb} G_i^b G_j^a = 0 \thinspace , \end{split} \end{equation}\] which reduces to \[\begin{equation} \label{eq:pse_ap1rog_real_orbitals} \begin{split} & g_{ai} \qty(1 - \qty(G^a_i)^2) + \qty[ 2 (h_{aa} - h_{ii}) + (g_{aa} - g_{ii}) + 2 \sum_{j \neq i}^{N_P} ( \gamma_{aj} - \gamma_{ij} ) ] G^a_i \\ & + \sum_{j \neq i}^{N_P} \qty(g_{ji} - g_{ja} G^a_i) G^a_j + \sum_{b = N_P + 1, b \neq a}^K \qty(g_{ab} - g_{ib} G^a_i) G^b_i + \sum_{b = N_P + 1, b \neq a}^K \sum_{j \neq i}^{N_P} g_{jb} G^a_j G^b_i = 0 \end{split} \end{equation}\] for real orbitals (De Baerdemacker 2017). After the \(N_P(K-N_P)\) equations (\(\eqref{eq:pse_ap1rog}\) or \(\eqref{eq:pse_ap1rog_real_orbitals}\)) are solved for the \(N_P(K-N_P)\) coefficients \(G_i^a\), the energy is given by \[\begin{align} E &= \sum_j^{N_P} \qty( 2 h_{jj} + \sum_k^{N_P} \gamma_{kj} + \sum_{b = N_P+1}^K g_{jb} G^b_j ) \\ &= \sum_j^{N_P} \qty( 2 h_{jj} + \sum_k^{N_P} (2 g_{kkjj} - g_{kjjk}) + \sum_{b = N_P + 1}^K g_{jbjb} G^b_j ) \label{eq:AP1roG_energy} \thinspace . \end{align}\] It is interesting to see (and quite expected as well, since AP1roG is pCCd) that the first two terms are exactly the RHF energy.
In the case of a weak interaction limit (De Baerdemacker 2017), i.e. \[\begin{equation} g_{jk} G_l^m \ll h_{aa} - h_{ii} \thinspace , \end{equation}\] the geminal coefficients can be approximated by \[\begin{equation} G_i^a = - \frac{g_{ai}}{2(h_{aa} - h_{ii})} \end{equation}\] at first order. This is especially useful as a better initial guess for the Newton algorithm.
In order to set up a procedure to solve these equations, let us first collect the \((K-N_P) \thinspace N_P\) (unknown) coefficients that appear in the AP1roG pSE equations (cfr. equation \(\eqref{eq:pse_ap1rog_real_orbitals}\)) in a coefficient vector \[\begin{equation} \vb{g} = (G_1^{N_P+1}, \ldots, G_1^K, G_2^{N_P+1}, \ldots, G_2^K, \ldots, G_i^a, \ldots, G_{N_P}^{N_P+1}, \ldots, G_{N_P}^K) \thinspace , \end{equation}\] and label them by indices \(\mu = 1, \ldots, (K-N_P) \thinspace N_P\).
When working in a computer environment, in which we will let \[\begin{align} i = 0, \ldots, N_P-1 \\ a = N_P, \ldots, K-1 \thinspace , \end{align}\] the conversion from the coefficient matrix \(\vb{G}\) (with elements \(G^a_i\), see equations \(\eqref{eq:coefficients_ap1rog_matrix_1}\) and \(\eqref{eq:coefficients_ap1rog_matrix_2}\)) to the coefficient vector \(\vb{g}\) (with elements \(g_\mu\)) can be done using \[\begin{equation} \label{eq:a_i_to_mu} \mu = (K - N_P) \thinspace i + a \thinspace , \end{equation}\] with \[\begin{equation} \mu = 0, \ldots, (K-N_P) N_P - 1 \thinspace . \end{equation}\] To go from the vector index \(\mu\) to the matrix indices \(i\) and \(a\), we can use the following formulas: \[\begin{align} &i = \lfloor{ \frac{\mu}{K - N_P} \rfloor} \\ &a = \mu \bmod(K - N_P) + N_P \thinspace . % use \bmod because I don't want a leading whitespace before mod \end{align}\]
Introducing the coordinate functions \[\begin{equation} \begin{split} f_i^a &: \mathbb{R}^{(K-N_P) N_P} \rightarrow \mathbb{R} \\ &: \vb{g} \mapsto g_{ai} \qty(1 - \qty(G^a_i)^2) + \qty[ 2 \sum_{j \neq i}^{N_P} ( \gamma_{aj} - \gamma_{ij} ) + (2 h_{aa} - 2 h_{ii}) + (g_{aa} - g_{ii}) ] G^a_i \\ &\phantom{: \vb{c} \mapsto} + \sum_{b = N_P + 1, b \neq a}^K \qty(g_{ab} - g_{ib} G^a_i) G^b_i + \sum_{j \neq i}^{N_P} \qty(g_{ji} - g_{ja} G^a_i) G^a_j + \sum_{b = N_P + 1, b \neq a}^K \sum_{j \neq i}^{N_P} g_{jb} G^a_j G^b_i \thinspace , \end{split} \end{equation}\] and \[\begin{equation} \begin{split} \vb{F} &: \mathbb{R}^{(K-N_P) N_P} \rightarrow \mathbb{R}^{(K-N_P) N_P} \\ &: \vb{g} \mapsto f_i^a(\vb{g}) \thinspace , \end{split} \end{equation}\] which is \[\begin{equation} \vb{F}(\vb{g}) = \begin{pmatrix} f_1^{N_P}(\vb{g}) \\ \vdots \\ f_i^a(\vb{g}) \\ \vdots \\ f_{N_P}^K(\vb{g}) \end{pmatrix} \thinspace , \end{equation}\] we can write the pSEs as \[\begin{equation} \label{eq:pse_f(g)} \vb{F}(\vb{g}) = \vb{0} \thinspace . \end{equation}\]
In order to solve equation \(\eqref{eq:pse_f(g)}\), we can use Newton’s method with the Jacobian matrix evaluated at \(\vb{g}\): \[\begin{equation} \vb{J}(\vb{g}) = \begin{pmatrix} \dfrac{\partial f_1}{\partial g_1}(\vb{g}) & \cdots & \dfrac{\partial f_1}{\partial g_{(K-N_P)N_P}}(\vb{g}) \\ \vdots & \ddots & \vdots \\ \dfrac{\partial f_{(K-N_P)N_P}}{\partial g_1}(\vb{g}) & \cdots & \dfrac{\partial f_{(K-N_P)N_P}}{\partial g_{(K-N_P)N_P}} (\vb{g}) \end{pmatrix} \thinspace , \end{equation}\] where the indices on \(f\) and \(g\) are labelled by \(\mu=1, \ldots, (K-N_P)N_P\) (cfr. equation \(\eqref{eq:a_i_to_mu}\)). For our specific problem, we can calculate \[\begin{equation} \begin{split} \pdv{f_i^a(\vb{g})}{G_k^c} = &- 2 \delta_{ac} \delta_{ik} g_{ai} G^a_i + \delta_{ac} \delta_{ik} \qty[ 2 \sum_{j \neq i}^{N_P} ( \gamma_{aj} - \gamma_{ij} ) + (2 h_{aa} - 2 h_{ii}) + (g_{aa} - g_{ii}) ] \\ & + \delta_{ik} \qty[ (1-\delta_{ac})(g_{ac} - g_{ic}G^a_i) - \delta_{ac} \sum_{b = N_P + 1, b \neq a}^K g_{ib} G^b_i ] \\ & + \delta_{ac} \qty[ (1-\delta_{ik})(g_{ki} - g_{ka}G^a_i) - \delta_{ik} \sum_{j \neq i}^{N_P} g_{ja} G^a_j ] \\ & + \delta_{ac} (1 - \delta_{ik}) \sum_{b = N_P + 1, b \neq a}^K g_{kb} G^b_i + \delta_{ik} (1 - \delta_{ac}) \sum_{j \neq i}^{N_P} g_{jc} G^a_j \thinspace . \end{split} \end{equation}\] Specifically, \[\begin{align} &c \neq a \land k \neq i \qquad \Rightarrow \qquad \pdv{f_i^a(\vb{g})}{G_k^c} = 0 \\ &c = a \land i \neq k \qquad \Rightarrow \qquad \pdv{f_i^a(\vb{g})}{G_k^a} = g_{ki} - g_{ka} G^{a}_i + \sum_{b = N_P + 1, b \neq a}^K g_{kb} G_i^b \\ &c \neq a \land i = k \qquad \Rightarrow \qquad \pdv{f_i^a(\vb{g})}{G_i^c} = g_{ac} - g_{ic} G^{a}_i + \sum_{j \neq i}^{N_P} g_{jc} G_j^a \end{align}\] and if both pairs of indices are equal (i.e. \(a=c \land i=k\)): \[\begin{equation} \pdv{f_i^a(\vb{g})}{G_i^a} = -2 g_{ai} G^a_i + 2 \sum_{j \neq i}^{N_P}(\gamma_{aj} - \gamma_{ij}) + (2 h_{aa} - 2 h_{ii}) + (g_{aa} - g_{ii}) - \sum_{b = N_P + 1, b \neq a}^K g_{ib} G_i^b - \sum_{j \neq i}^{N_P} g_{ja} G_j^a \thinspace . \end{equation}\]
\
If the recomputation of the Jacobian at every Newton step should be too hard, we could also use Broyden’s method.
The AP1roG energy is \[\begin{equation} E = % \sum_j^{N_P} \qty( % 2 h_{jj} + \sum_k^{N_P} \gamma_{kj} % + \sum_{b = N_P + 1}^K g_{jb} G^b_j % ) % \end{equation}\] and the PSEs are (in the real case) \[\begin{equation} \begin{split} & g_{ai} \qty(1 - \qty(G^a_i)^2) % + \qty[ % 2 (h_{aa} - h_{ii}) % + (g_{aa} - g_{ii}) % + 2 \sum_{j \neq i}^{N_P} ( % \gamma_{aj} - \gamma_{ij} % ) % ] G^a_i \\ & + % \sum_{j \neq i}^{N_P} \qty( % g_{ji} - g_{ja} G^a_i % ) G^a_j % + \sum_{b = N_P + 1, b \neq a}^K \qty( % g_{ab} - g_{ib} G^a_i % ) G^b_i % + \sum_{b = N_P + 1, b \neq a}^K \sum_{j \neq i}^{N_P} % g_{jb} G^a_j G^b_i % = 0 % \thinspace . \end{split} \end{equation}\]