Numerical methods for computing variationally optimized molecular orbitals within the HartreeFock approximation are augmented to include correlation functionals of the density in the energy and the numerical methods for carrying this out are described. The approach is applied explicitly to the ColleSalvetti correlation energy functional. It is found that the gradient terms in the ColleSalvetti functional present numerical problems associated with the low-density behavior, but also that they make a relatively small contribution to the correlation energy. In the three cases considered, HF, H2O and N2, it is found that the ColleSalvetti correction considerably underestimates the correlation energies obtained in coupled-cluster theory.