We present a new algorithm for performing ab initio solution phase geometry optimizations. The procedure is based on the self consistent-reaction-field method developed in our laboratory which combines electronic structure calculations with a finite element formulation of the continuum electrostatics problem. A gradient for the total solution phase free energy is obtained by combining different contributions from the gradient of the classical polarization free energy and the derivatives of the quantum mechanical energy. The method used in obtaining the classical gradient is based on exact linear algebra relations and a Green function formalism due to Handy and Schaefer. Both the classical and quantum mechanical gradients are validated by comparison with energy finite differences. The result of applications to a number of small organic compounds are discussed. Comparisons between the predicted location and depth of the various solution phase minima of the Ramachandran map for the alanine dipeptide and those reported by Gould et al. are also presented.