Variational optimization of the second order density matrix derived from a seniority-zero wave function

We perform a direct variational determination of the second-order (two-particle) density matrix corresponding to a many-electron system, under a restricted set of the two-index $N$-representability P-, Q- and G-conditions. In addition, we impose the constraints that the two-particle density matrix must be derivable from a doubly-occupied many-electron wave function, i.e.\ a singlet wave function for which the Slater determinant decomposition only contains determinants in which spatial orbitals are doubly occupied. We rederive the two-index conditions first found in [F. Weinhold and E. Wilson J. Chem. Phys. 46, 2752 (1967)] and [F. Weinhold and E. Wilson, J. Chem. Phys. 47, 2298 (1967)] and apply them to various systems (hydrogen chains, $\text{He}$, $\text{H}_2$, dissociation of diatomic molecules). Compared to the general case, the structure of doubly-occupied two-particle density matrices causes the associate semidefinite program to have a very favorable scaling as $L^3$, where $L$ is the number of spatial orbitals. Since the doubly-occupied Hilbert space depends on the choice of the orbitals, variational calculation steps of the two-particle density matrix are interspersed with orbital-optimization steps (based on Jacobi rotations in the space of the spatial orbitals). We also point to the importance of symmetry breaking of the orbitals when performing calculations in a doubly-occupied framework.