We present full variational optimization of the natural orbitals and their occupation numbers for two functionals of the one-electron reduced density matrix along the entire dissociation coordinate of the H-2, Li-2 and LiH molecules. The corrected Hartree functional is characterized by an exchange-correlation part in which the square roots of the NO occupation numbers enter. It performs well in atomic and molecular systems, when near-exact Nos and NO occupation numbers are used. However, it is here shown to fail the test of variational stability, yielding too low energies at equilibrium geometry and in particular at long bond lengths. The corrected Hartree-Fock functional on the other hand improves little or not at all over Hartree-Fock around the equilibrium bond length. The search for a satisfactory expression for the dependence of the natural orbital functionals on the Nos and NO occupation numbers has to be continued. (C) 2002 Elsevier Science B.V. All rights reserved.