To date, conventional force fields such as GAFF, UFF, CHARMM, and OPLS-AA have been the main workhorses in materials science and biological simulations. All these force fields partition the total energy into inter-and intra-molecular parts based on the bonding topology of the molecule. Intermolecular interactions are described by the atomic point charges and simple pairwise-additive van der Waals (VDW) potentials. The intramolecular interactions are decomposed into a direct sum of contributions from all internal coordinates (i.e., bond lengths, angles, dihedrals etc.). These force fields can be computed in fast speed, making them primary models to use in industry to study large molecular systems. However, the mathematical form of the conventional force fields is heavily approximated, neglecting many important physics including polarization, charge penetration, and many-body dispersion etc. Therefore, conventional force fields are hard to generalize, and are essentially effective models that need to be tailored for specific systems or applications. The force field parameters are often obtained through a top-down approach, fitted to reproduce the experimental macroscopic properties of interests. Such parameterization process relies heavily on human intervention, thus cannot be done in high throughput. Consequently, force field parameterization is often the bottleneck of research when studying new systems with a large number of parameters. There are efforts trying to automate this process, such as ForceBalance and OpenFF, but the differentiation of macroscopic properties with respect to force field parameters is a nontrivial task. In these works, such differentiation is still conducted using the finite difference method, resulting in limited efficiency and poor scalability. Considering how widely conventional force fields are used in modern industry, their automatic optimization remains an important challenge to be solved.