In this study, an electrostatic potential (ESP) fitting method using constrained spatial electron density (cSED) expanded with preorthogonal natural atomic orbitals (pNAOs) was proposed. In this method, the electron density of a molecule is divided into spherical atom-centered electron densities and the expansion coefficient is determined to reproduce the ESP around the molecule. Our method was then applied to two systems: (i) a hydration reaction of cis-platin and (ii) a variety of organic/inorganic molecules. By evaluating the atomic charges along the hydration reaction, our method showed good conformational transferability, which cannot be obtained using conventional ESP fitting methods. Moreover, we successfully obtained the hydration structure along the reaction by coupling our method with a reference interaction site model (RISM). Reasonable data were obtained not only for organic molecules but also for inorganic molecules. This success came from the introduction of pNAOs as auxiliary basis sets in the charge fitting.