This article derives theoretical foundations of force field functional theory (FFFT). FFFT studies topics related to the functional representation of nonreactive forcefields to achieve various desirable properties such as: (a) formal exactness of the forcefield's energy functional under certain conditions, (b) a formally exact ansatz separating the bonded potential energy from the nonbonded potential energy within a bonded cluster in a way that enables bonded parameters to be optimized using linear regression instead of requiring nonlinear regression, (c) the potential energy's continuous differentiability to various orders with respect to energetically accessible internal coordinate displacements within a subdomain defined by one electronic ground state, (d) forcefield design that guarantees the reference ground-state geometry is exactly reproduced as an equilibrium structure on the forcefield's potential energy landscape, (e) reasonably accurate and broadly applicable frugal model potentials, (f) computationally efficient embedded feature selection that identifies and removes unimportant forcefield terms, (g) well-designed methods to parameterize the forcefield from quantum-mechanically-computed and (optionally) experimental reference data, and (h) forcefields that approximately reproduce experimentally-measured properties. This article also introduces: (1) an angle-bending model potential that more accurately describes physical dynamics and is continuously differentiable to all orders with respect to internal coordinate displacements even when the bond angle is linear (i.e., θ = π (180°)) and (2) a first-principles-derived stretch potential that accurately describes short-range Pauli repulsion and the long-range bond dissociation energy. This new angle-bending potential gave good agreement to CCSD quantum-chemistry calculations for CaH2, CO2, H2O, HNO, Li2O, NO2, NS2, SF2, SiH2, and SO2 molecules. This new bond-stretch potential reproduced the first 12+ and 30+ vibrational energy levels of H2 and O2 molecules, respectively, within a few percent of experimental values. Studying the C-F bond stretch in C6F6 as an example, the new ansatz (item (b) above) reduced sensitivity of the optimized force constant's value to choice of nonbonded interaction parameters by an order of magnitude compared to the old ansatz. Normal mode analysis of optimized flexibility models for CO2, H2O, HNO, and SO2 molecules yielded vibrational transition frequencies within a few percent of experimental values. These results demonstrate advantages of this new approach.