We present a new method for generating global or semiglobal potential energy surfaces in the presence of an electrostatic potential; the new method can be used to model chemical reactions in solution or in an enzyme, nanocavity, or other chemical environment. The method extends the multiconfiguration molecular mechanics method so that the energy depends on the electrostatic potential at each atomic center. The charge distribution of the system can also be calculated. We illustrate the method by applying it to the symmetric bimolecular reaction Cl- + CH3Cl → ClCH3 + Cl- in aqueous solution, where the potential energy information is obtained by the combined density functional and molecular mechanical method, that is, by the combined quantum mechanical and molecular mechanical method (QM/MM) with the QM level being density functional theory. It is found that we can describe a semiglobal potential energy surface in aqueous solution with electronic structure information obtained entirely in the gas phase, including the linear and quadratic responses to variations in the electrostatic potential distribution. The semiglobal potential energy surface calculated by the present method is in good agreement with that calculated directly without any fitting.