A new method for the treatment of many-body polarization effects in fluid systems is presented, making use of hybrid QM/MM techniques and the semiempirical wave function. In this approach, the electronic structure of each solvent molecule is represented by an antisymmetric determinant wave function, which can be determined by fully converged variational QM calculations at each step of the classical trajectory or can be treated as dynamic variables along with the nuclear coordinates in molecular dynamics simulations. In determining the molecular wave functions, a hybrid QM/MM method is used. This molecular orbital derived empirical potential for liquid simulations (MODEL) is analogous to the fluctuating charge model introduced by Berne and co-workers; however, the MODEL potential is a quantum chemical model as opposed to a purely empirical charge equilibration scheme employed in the FC model. The method may also be extended for the treatment of large molecular systems. The capability of the MODEL potential for describing bimolecular interactions was illustrated for bimolecular complexes. Reasonable agreement with ab initio results suggests that the MODEL potential may be parametrized for liquid simulations.