The methodology for treatment of proton transfer processes by density matrix evolution (DME) with inclusion of many excited states is presented. The DME method (Berendsen, H. J. C.; Mavri, J. J. Phys. Chem. 1993, 97, 13464) that simulates the dynamics of quantum systems embedded in a classical environment is extended to nonorthogonal basis sets. The method is applied to calculations of the rate of proton transfer in the double-well intramolecular hydrogen bond in the hydrogen malonate ion in aqueous solution. Five optimized basis functions were applied, giving rise to five vibrational levels of the proton. A nearly 30-fold increase in the proton transfer rate relative to the earlier treatment using only two basis functions was observed, indicating the importance of flexible, optimized basis sets in the proton transfer calculations. Possible applications of the method to biologically interesting processes such as proton transfer in enzymatic reactions are discussed.