The combination of all 4 models comprise a unique range of approximations used in molec- energy plots, as it has neither an explicit form of polarization nor a different frozen energy E frz from it's total interaction energy E int . Multiplet results are not shown here but are present in the supporting material as their plots are almost identical to the H-bonding results in Figure . All 4 models have comparable accuracy when computing the total interaction energy E int between nucleobases, with Drude-2017 having the largest median (-9.16kJ/mol) from zero for it's H-bonded configurations, but recovering some accuracy for stacked geometries where all models has less than a 1 kCal/mol median difference. While the spread of all distributions shown here have magnitudes close to their median values, it does improve when only looking at structures taken from the first binned category (See Supporting Info). The majority of the outliers are therefore mostly from geometries that are less likely to be found in NMR resolved structures.