Real singlet orbital optimization in the projected Schrödinger equation framework
The orbital rotation generator for real singlet rotations has the following form: \[\begin{equation} \hat{\kappa} = \sum_{p>q}^K \kappa_{pq} \hat{E}^-_{pq} \thinspace . \end{equation}\] Using this form, we can write the variational conditions on the orbital optimization Lagrangian as (\(p>q, \forall \boldsymbol{\eta}\)): \[\begin{align} \require{physics} & \eval{ \pdv{ \mathscr{L}(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}_0, \boldsymbol{\lambda}) }{\lambda_a} }_{\boldsymbol{\lambda}^\star} = 0 = f_a(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}_0) \\ % & \eval{ \pdv{ \mathscr{L}(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}, \boldsymbol{\lambda}^\star) }{\kappa_{pq}} }_{\boldsymbol{\kappa}_0} = 0 \\ % & \eval{ \pdv{ \mathscr{L}(\boldsymbol{\eta}, \vb{p}, \boldsymbol{\kappa}_0, \boldsymbol{\lambda}^\star) }{p_i} }_{\vb{p}^\star} = 0 = \eval{ \pdv{ E(\boldsymbol{\eta}, \vb{p}, \boldsymbol{\kappa}_0) }{p_i} }_{\vb{p}^\star} + \sum_a^S \lambda_a^\star \qty( \eval{ \pdv{ f_a(\boldsymbol{\eta}, \vb{p}, \boldsymbol{\kappa}_0) }{p_i} }_{\vb{p}^\star} ) \thinspace . \end{align}\] Note that we are employing standard orbital optimization techniques: evaluating the orbital rotation parameters \(\boldsymbol{\kappa}\) at \(\boldsymbol{\kappa}_0 = \vb{0}\) refers to the situation that we employ the current orbitals at every iteration step of the orbital optimization.
The stationary condition with respect to the orbital rotation parameters \(\kappa_{pq}\) means that the orbital gradient of the Lagrangian should vanish, which is commonly referred to as the generalized Brillouin theorem. Let us now try to find an expression for the orbital gradient of the Lagrangian. By using the BCH expansion, we can find the orbital gradient of the energy function: \[\begin{equation} \eval{ \pdv{ E(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}) }{\kappa_{pq}} }_{\boldsymbol{\kappa}_0} = \frac{ \matrixel{0}{ \comm*{ \hat{E}^-_{pq} }{ \hat{\mathcal{H}}(\boldsymbol{\eta}) } }{\Psi(\vb{p}^\star)} }{ \braket{0}{\Psi(\vb{p}^\star)} } \end{equation}\] and the orbital gradient of the PSE functions: \[\begin{equation} \eval{ \pdv{ f_a(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}) }{\kappa_{pq}} }_{\boldsymbol{\kappa}_0} = \matrixel{a}{ \comm*{ \hat{E}^-_{pq} }{ \hat{\mathcal{H}}(\boldsymbol{\eta}) } }{\Psi(\vb{p}^\star)} \\ - \frac{ \braket{a}{\Psi(\vb{p}^\star)} }{ \braket{0}{\Psi(\vb{p}^\star)} } \matrixel{0}{ \comm*{ \hat{E}^-_{pq} }{ \hat{\mathcal{H}}(\boldsymbol{\eta}) } }{\Psi(\vb{p}^\star)} \thinspace . \end{equation}\] Then, by adequately defining the response Fockian matrix in terms of the response density matrices (in the real case): \[\begin{equation} 2 \thinspace \mathscr{F}^{\text{r}}_{pq} ( \hat{\mathcal{H}}(\boldsymbol{\eta}) ) = \sum_r^K h_{qr}(\boldsymbol{\eta}) \qty( D^{\text{r}}_{pr} + D^{\text{r}}_{rp} ) + \sum_{rst}^K g_{qrst} (\boldsymbol{\eta}) \qty( d^{\text{r}}_{prst} + d^{\text{r}}_{rpst} ) \thinspace , \end{equation}\] we can calculate the orbital gradient as: \[\begin{equation} \eval{ \pdv{ \mathscr{L}(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}, \boldsymbol{\lambda}^\star) }{\kappa_{pq}} }_{\boldsymbol{\kappa}_0} = 2 \qty[ \mathscr{F}^{\text{r}}_{pq} ( \hat{\mathcal{H}}(\boldsymbol{\eta}) ) - \mathscr{F}^{\text{r}}_{qp} ( \hat{\mathcal{H}}(\boldsymbol{\eta}) ) ] \thinspace , \end{equation}\] which has the same form as regular variational approaches.
The orbital Hessian then becomes (with \(p>q, r>s\)): \[\begin{equation} \eval{ \pdv{ \mathscr{L}(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}, \boldsymbol{\lambda}^\star) }{\kappa_{pq}}{\kappa_{rs}} }_{\boldsymbol{\kappa}_0} = (1 + P_{pr})(1 - P_{pq})(1 - P_{rs}) \mathscr{G}^{\text{r}}_{pqrs} ( \hat{\mathcal{H}}(\boldsymbol{\eta}) ) \thinspace , \end{equation}\] in which we have adequately defined the response super-Fockian matrix: \[\begin{align} 2 \thinspace \mathscr{G}^{\text{r}}_{pqrs} ( \hat{\mathcal{H}}(\boldsymbol{\eta}) ) = &2 \delta_{qr} \mathscr{F}^{\text{r}}_{ps} ( \hat{\mathcal{H}}(\boldsymbol{\eta}) ) - h_{qr} ( D^{\text{r}}_{ps} + D^{\text{r}}_{sp} ) \\ &+ \sum_{tu}^K \Big( g_{qtsu} ( d^{\text{r}}_{ptru} + d^{\text{r}}_{tpur} ) - g_{qrtu} ( d^{\text{r}}_{pstu} + d^{\text{r}}_{sput} ) - g_{qtur} ( d^{\text{r}}_{ptus} + d^{\text{r}}_{tpsu} ) \Big) \thinspace . \end{align}\]
By virtue of the evaluation at \(\boldsymbol{\kappa}_0\), which means that the exponential vanishes, the first two conditions on the orbital optimization Lagrangian are equal to their non-orbital optimized counterparts. We can therefore propose the following iterative and sequential orbital optimization algorithm (in the current orbital basis \(\boldsymbol{\kappa}_0\)):
solve the PSEs for the wave function parameters \(\vb{p}^\star\);
subsequently solve the linear equations for the Lagrange multipliers \(\boldsymbol{\lambda}^\star\), from which we can calculate the response-1- and 2-DMs;
use the response-1- and 2-DMs to check if the current Lagrangian orbital gradient is zero. If it is not, we could for example use a Newton step \[\begin{equation} \eval{ \pdv[2]{ \mathscr{L}(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}, \boldsymbol{\lambda}^\star) }{\boldsymbol{\kappa}} }_{\boldsymbol{\kappa}_0} \boldsymbol{\kappa} = - \eval{ \pdv{ \mathscr{L}(\boldsymbol{\eta}, \vb{p}^\star, \boldsymbol{\kappa}, \boldsymbol{\lambda}^\star) }{\boldsymbol{\kappa}} }_{\boldsymbol{\kappa}_0} \end{equation}\] and use the parameter update \(\boldsymbol{\kappa}\) to rotate the orbital basis into a better one and start a next iteration.
For the specific case of coupled cluster theory, we can find similar formulas for the orbital optimization coupled-cluster Lagrangian by Helgaker (T. Helgaker, Jørgensen, and Olsen 2000). Bozkaya (Bozkaya et al. 2011), has also formulated a (quadratically convergent) orbital optimization algorithm for CCD, but we should note that their definitions for the response density matrices and response Fockian matrix are slightly different.