The optimization of carbonaceous environmental sensors is currently limited by empirical statistical models that fail at quantum-mechanical boundaries. In this study, the standard Response Surface Methodology (RSM) diverges significantly in highly non-linear van der Waals regimes (MAE > 9 eV); to address this, a Density Functional Theory with D3 dispersion correction (DFT-D3) and an explainable-AI framework are employed. By computing a deterministic 500-point 4D phase space ( \(\:X,\:Y,\:Z,\:\theta\:\) ) for caffeine physisorption on nanobiochar, a Gaussian Process Regression surrogate with a Matern-5/2 kernel (cross-validated MAE = 0.006 eV, below chemical accuracy) was trained to efficiently identify the lowest-energy configurations within the sampled phase space. Symbolic Regression (PySR) was subsequently applied to derive a continuous analytical equation that explicitly couples vertical adsorbate displacement to the lateral periodic registry of the \(\:s{p}^{2}\) lattice. Electronic structure analysis of the lowest-energy configuration confirms an eclipsed \(\:\pi\:-\pi\:\) stacking architecture at an equilibrium distance of 3.4 Å. Out-of-plane \(\:{p}_{z}\) orbital hybridization induces localized donor-acceptor charge transfer, resulting in a 1.06 eV shift in the contact potential difference (increasing the local work function from 4.14 eV to 5.20 eV). The absolute magnitude of this intrinsic electronic binding enthalpy exceeds the estimated classical entropic penalty of surface confinement, indicating thermodynamically favorable adsorption ( \(\:\varDelta\:{G}_{ads}<0\) ) at 298.15 K. This methodology provides a physically interpretable blueprint for predicting the selectivity and transduction mechanisms of molecular adsorbates on carbonaceous substrates.
