A variety of classical (pairwise) force fields have been utilized to elucidate the underlying mechanism but failed to reproduce the frequency-dependent dielectric responses quantitatively. For instance, nonpolarizable models (e.g., SPC/E and TIP4P/2005 ) cannot reproduce the static dielectric constant and the cooperative relaxation time. Including the bond flexibility in the pairwise model (e.g., SPC/Fw and q-TIP4P/F ) improves predicting the static dielectric constant. However, these models lack many-body interactions ; the three-body interaction energy contributes to approximately 15 -20% of the total interaction energy . While MD simulations with pair potentials can yield qualitatively acceptable results, obtaining "the right results for the right reasons" by reflecting the many-body interactions is essential for interpreting the underlying molecular mechanisms .