A quantum Monte Carlo method of determining Jastrow–Slater and correlated multideterminant wave functions for which the energy is stationary with respect to variations in the single-particle orbitals is presented. A potential is determined by a least-squares fitting of fluctuations in the energy with a linear combination of one-body operators. This potential is used in a self-consistent scheme for the orbitals whose solution ensures that the energy of the correlated wave function is stationary with respect to variations in the orbitals. The method is feasible for atoms, molecules, and solids and is demonstrated for the beryllium, carbon, and neon atoms and for the solid diamond.