We present an extended methodology of the quantum mechanical/molecular mechanical framework with Metropolis Monte Carlo method (QM/MM-MC), in which close surrounding solvent molecules are included in QM subsystem to evaluate quantitatively the solvation free energy change. The neutral (N) and zwitterionic (Z) structures of glycine in aqueous solution were optimized in the combination procedure of simulated annealing with MC at MM level and QM/MM-vib geometry optimization. Helmholtz energy change between N and Z forms of glycine in aqueous solution was calculated stepwise with expanded QM/MM-MC method, and the computational result was in good agreement with the experimental value.