In the past few decades, molecular mechanics has become an important tool for interpreting and predicting properties of materials and biomolecules. Under Born-Oppenheimer approximation, all molecular mechanics simulations rely on a potential energy surfaces (PES) to describe the electronic interactions between atoms. The quality of the PES directly controls the reliability of the simulation, and it can be either computed on-the-fly using ab initio methods or approximated by a simplified "force field" model. Either way, it remains a great challenge to obtain a reliable PES due to its high dimensionality. High-level correlated wavefunction (CW) methods beyond density functional theory (DFT) are haunted by the "curse of dimensionality", which makes them prohibitively expensive for systems more than a few tens of atoms. It is also nontrivial to construct a robust force field model, as it requires a deep understanding to the physical nature of all the interactions in the system. For molecular systems with weak nonbonding interactions, the nonempirical force field development is especially difficult because coupled-cluster level of theory is often needed to provide enough accuracy. The computational cost of CCSD(T) scales as 7 () N  , so the reference data needed for the force field fitting is very hard (if not impossible) to obtain for large systems.