Edmiston-Ruedenberg localization
This scheme (Weinstein, Pauncz, and Cohen 1971) is inspired by the fact that the two-electron interaction energy of a single Slater determinant, \[\begin{equation} \require{physics} E_\text{int} % = \sum_{ij}^{N_P} (2 J_{ij} - K_{ij}) % = \sum_{ij}^{N_P} (2 g_{iijj} - g_{ijji}) \end{equation}\] is invariant with respect to unitary transformations of the occupied orbitals, but the term \[\begin{equation} \label{eq:Edmiston-Ruedenberg_localization_index} D = \sum_i^{N_P} J_{ii} = \sum_i^{N_P} g_{iiii} \end{equation}\] is not. This term seems to be related to a self-interaction energy of an electron occupying the orbital \(\phi_i\), and maximizing its value can therefore be used as a criterion for producing localized (occupied) orbitals from the canonical RHF orbitals and is coined the sum of orbital self-interaction energies.
The full transformation is then written as \[\begin{align} D' &= \sum_i^{N_P} g'_{iiii} \\ &= \sum_{ijklm}^{N_P} % U^*_{ji} U_{ki} U^*_{li} U_{mi} % g_{jklm} % \thinspace , \end{align}\] or we can make the transformation-dependence clear by writing \[\begin{equation} D(\boldsymbol{\kappa}) % = \sum_{ijklm}^{N_P} % \qty[ \exp(-\boldsymbol{\kappa}) ]^*_{ji} % \qty[ \exp(-\boldsymbol{\kappa}) ]_{ki} % \qty[ \exp(-\boldsymbol{\kappa}) ]^*_{li} % \qty[ \exp(-\boldsymbol{\kappa}) ]_{mi} % g_{jklm} % \thinspace , \end{equation}\] in which the orbital rotation parameters are those appearing in the generator for real, singlet occupied-occupied orbital rotations \[\begin{equation} \hat{\kappa} = % \sum_{i>j}^{N_P} \kappa_{ij} \hat{E}^-_{ij} % \thinspace . \end{equation}\]
Using Jacobi rotations, we can write \[\begin{equation} D' = D + A + B \cos(4 \theta) + C \sin(4 \theta) \thinspace , \end{equation}\] in which \[\begin{align} & A = \frac{1}{4} ( % 2 g_{iijj} + 4 g_{ijij} - g_{iiii} - g_{jjjj} % ) \\ & B = \frac{1}{4} ( % g_{iiii} + g_{jjjj} - 2 g_{iijj} - 4 g_{ijij} % ) \\ & C = g_{ijjj} - g_{iiij} % \thinspace , \end{align}\] in which \(i>j\) are two orbital indices of orbitals that are occupied in the Hartree-Fock Slater determinant. Its maximum is given by \[\begin{equation} D_\text{max}' = D + A + \sqrt{B^2 + C^2} \thinspace , \end{equation}\] for the maximizing angle \[\begin{equation} \theta_\text{max} = % \frac{1}{4} \arctan \qty( % \frac{B}{\sqrt{B^2 + C^2}}, % \frac{C}{\sqrt{B^2 + C^2}} % ) % \thinspace . \end{equation}\] We can then formulate an algorithm to localize a given spatial orbital basis using the Jacobi-Edmiston-Ruedenberg localization method. We also require the number of electron pairs \(N_P\) and a threshold \(\epsilon\) to check for convergence.As opposed to using Jacobi rotations to maximize the Edmiston-Ruedenberg localization index, we can calculate the derivatives with respect to the parameters \(\boldsymbol{\kappa}\). For the gradient, we can write \[\begin{equation} \vb{D}^{(1)}_{ij} = % \eval{ % \pdv{ % D(\boldsymbol{\kappa}) % }{ \kappa_{ij} } % }_{ \boldsymbol{\kappa} = \vb{0} } % = \sum_{cd}^K % \eval{ % \pdv{ D(\vb{U}) }{ U_{cd} } % }_{ \vb{U} = \vb{I} } % \qty[ % \eval{ % \pdv{ \vb{U} }{ \kappa_{ij} } % }_{ \boldsymbol{\kappa} = \vb{0} } % ]_{cd} \end{equation}\] and for the Hessian: \[\begin{equation} \begin{split} \vb{D}^{(2)}_{ijkl} % = \eval{ % \pdv{ % D(\boldsymbol{\kappa}) % }{ \kappa_{ij} }{ \kappa_{kl} } % }_{ \boldsymbol{\kappa} = \vb{0} } % = & \sum_{cd}^K % \eval{ % \pdv{ D(\vb{U}) }{ U_{cd} } % }_{ \vb{U} = \vb{I} } % \qty[ % \eval{ % \pdv{ \vb{U} }{ \kappa_{ij} }{ \kappa_{kl} } % }_{ \boldsymbol{\kappa} = \vb{0} } % ]_{cd} \\ &+ \sum_{abcd}^K % \eval{ % \pdv{ D(\vb{U}) }{ U_{ab} }{ U_{cd} } % }_{ \vb{U} = \vb{I} } % \qty[ % \eval{ % \pdv{ \vb{U} }{ \kappa_{ij} } % }_{ \boldsymbol{\kappa} = \vb{0} } % ]_{ab} % \qty[ % \eval{ % \pdv{ \vb{U} }{ \kappa_{kl} } % }_{ \boldsymbol{\kappa} = \vb{0} } % ]_{cd} % \thinspace . \end{split} \end{equation}\] For the first derivative, we can calculate: \[\begin{equation} \pdv{ D(\vb{U}) }{ U_{cd} } % = 4 \sum_{ijk}^{N_P} % U_{id} U_{jd} U_{kd} % g_{cijk} % \thinspace , \end{equation}\] which becomes (evaluated at \(\vb{U} = \vb{I}\)): \[\begin{equation} \eval{ % \pdv{ D(\vb{U}) }{ U_{cd} } % }_{ \vb{U} = \vb{I} } % = 4 g_{cddd} % \thinspace , \end{equation}\] such that the gradient has a pretty simple expression: \[\begin{equation} \vb{D}^{(1)}_{ij} % = \eval{ % \pdv{ D(\boldsymbol{\kappa}) }{ \kappa_{ij} } % }_{ \boldsymbol{\kappa} = \vb{0} } % = 4 (g_{jiii} - g_{ijjj}) % \thinspace . \end{equation}\] For the second derivative, we have \[\begin{equation} \pdv{ D(\vb{U}) }{ U_{ab} }{ U_{cd} } % = 4 \delta_{bd} \sum_{ij}^K % U_{ib} U_{jb} (2 g_{ciaj} + g_{caij}) % \thinspace , \end{equation}\] which becomes (evaluated at \(\vb{U} = \vb{I}\)): \[\begin{equation} \eval{ % \pdv{ D(\vb{U}) } {U_{ab} }{ U_{cd} } % }_{ \vb{U} = \vb{I} } % = 4 \delta_{bd} (2 g_{cbab} + g_{cabb}) % \thinspace , \end{equation}\] such that we get for the Hessian: \[\begin{equation} \begin{split} \vb{D}^{(2)}_{ijkl} = % & \delta_{ik} ( % -2 g_{jlll} - 2 g_{ljjj} + 8 g_{liji} + 4 g_{ljii} % ) \\ &+ \delta_{jk} ( % 2 g_{illl} + 2 g_{liii} - 8 g_{ljij} - 4 g_{lijj} % ) \\ &+ \delta_{il} ( % 2 g_{jkkk} + 2 g_{kjjj} - 8 g_{kiji} - 4 g_{kjii} % ) \\ &+ \delta_{jl} ( % -2 g_{ikkk} - 2 g_{kiii} + 8 g_{kjij} + 4 g_{kijj} % ) % \thinspace . \end{split} \end{equation}\]