We proposed a new fitting approach of electrostatic potential (ESP) for stable quantum mechanical (QM) calculations using the reference interaction site model. In this approach, the approximated density matrix is fitted so that the ESP computed from QM calculations can be reproduced. We introduced two conditions: conservation of the number of electrons and positive semi-definite condition of the electron density matrix. This approach was introduced in a linear response approximation of the reference interaction site model self-consistent field explicitly including the spatial electron density distribution. By considering the two conditions, we overcame the instability inherent in a previous approach. (C) 2013 Elsevier B. V. All rights reserved.