  Ground state eigenvectors of the reduced Bardeen-Cooper-Schrieffer
Hamiltonian are employed as a wavefunction ansatz to model strong electron
correlation in quantum chemistry. This wavefunction is a product of
weakly-interacting pairs of electrons. While other geminal wavefunctions may
only be employed in a projected Schr\"{o}dinger equation, the present approach
may be solved variationally with polynomial cost. The resulting wavefunctions
are used to compute expectation values of Coulomb Hamiltionans and we present
results for atoms and dissociation curves which are in agreement with
doubly-occupied configuration interaction (DOCI) data. The present approach
will serve as the starting point for a many-body theory of pairs, much as
Hartree-Fock is the starting point for weakly-correlated electrons.
